This script is used to run LUKE for the runaway problem by J. Decker, Y. Peysson and E. Nilsson
0001 % 0002 % This script is used to run LUKE for the runaway problem 0003 % 0004 % by J. Decker, Y. Peysson and E. Nilsson 0005 % 0006 clear all 0007 clear mex 0008 clear functions 0009 close all 0010 warning('off') 0011 % 0012 % *********************** This part must be specified by the user ***************************** 0013 % 0014 % printing option 0015 % 0016 permission = test_permissions_yp; 0017 % 0018 if ~permission 0019 disp('Please move the script to a local folder where you have write permission before to run it') 0020 return; 0021 end 0022 % 0023 id_simul = 'Runaway_alpha3_nosync';%_coll2%Simulation ID 0024 % id_simul = 'Runaway_alpha3_sync_0p6';%_coll2_norm1_Zi4%Simulation ID 0025 p_opt = -1;% printing and saving options : (-1) do nothing (0) print figures (1) print figures and save figures and results (2) save figures and results 0026 % 0027 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants 0028 % 0029 % calculation mode : (0) quasi steady state (1) time evolution 0030 % 0031 opt.timevol = 1; 0032 % 0033 % normalized electric field (with respect to the Critical field Ec) 0034 % 0035 alpha = 3;%3.25;% E/Ec 0036 % 0037 % temperature betath2 = Te/mc2 0038 % (betath2 = 10^(-6) is validated for NR limit) 0039 % 0040 betath2 = 0.01;% Te/mc2 0041 % 0042 % the ratio of the syncrotron reaction force to the collisional drag scales 0043 % as taurinv = 1/taur 0044 % -> set taurinv = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0045 % 0046 % note : 0047 % - without the syncrotron reaction force, and with a single cold ion 0048 % species of charge Zi, the normalized runaway rate - with respect to the 0049 % electron density and ee collision time - only depends upon Zi, alpha and betath2 0050 % 0051 % - with the syncrotron reaction force, an additional parameter is 0052 % provived by taurinv 0053 % 0054 % - the link between wpr=wc2/wp2 and taurinv is : 0055 % taurinv = wpr*(2/3)/ln(Lambda) 0056 % 0057 % - the magnetic field being fixed, the density is adjusted in LUKE according to bhat 0058 % 0059 taurinv = 0.6;%0.3;%1;%0.7523;%0.9;% 0060 % 0061 opt.synchro_mode = 0; 0062 % opt.synchro_mode = 1; 0063 % 0064 % ion charge (single species, cold ions) 0065 % 0066 Zi = 1;%4;% 0067 % 0068 Emin = NaN;%runaway energy threshold in MeV (limit where electrons are "counted") 0069 % Emax = (sqrt(20^2*betath2+1)-1)*mc2/1e3;%1;%grid energy threshold in MeV [(sqrt(pnmax_S^2*betath2+1)-1)*mc2/1e3 = 0.63163 corresponds to pnmax_S = 20 for betath = 0.1] 0070 % 0071 % collision mode : 0072 % (0) : Relativistic Maxwellian background 0073 % (1) : High-velocity limit 0074 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0075 % 0076 opt.coll_mode = 0;%2;% 0077 opt.bounce_mode = 0;%1;% bounce-averaged calculation 0078 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0079 % 0080 % normalization mode : 0081 % (0) : "free" mode 0082 % (1) : normalization of f0 maintained to its initial (maxwellian) value. DO NOT use in the runaway problem 0083 % (2) : normalization of f0 maintained to its value at the previous time. Use only in the runaway problem with (+1i) 0084 % (3) : injection at center of particles lost at boundary (use only for stream function) 0085 % (... + 1i) : compensate for numerical error after each time step, acounting for boundary losses. 0086 opt.norm_mode_f = 3;%1i;%1;%0;%3+ 0087 % 0088 % time grid (normalized to collision time) -> you can specify : 0089 % - linear arbitrary grid with tn array : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0 0090 % - linear arbitrary grid with dtn array : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10); 0091 % - linear grid with tnmax and grid step : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000; 0092 % - linear grid with tnmax and number of grid steps : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i; 0093 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10; 0094 % 0095 % opt.tnmin = 0;% > 0 for log scale 0096 % opt.tn = 10000;% 10000 is time for asymptotic solution 0097 % opt.dtn = 1000;% 10 time steps required for accurate runaway solution - see rundke_dtn 0098 % graph.itdisp = [1,5,10];% display times - only used if opt.timevol > 0 0099 % 0100 % opt.tnmin = 1;% > 0 for log scale 0101 % opt.tn = 1e5;% 10000 is time for asymptotic solution 0102 % opt.dtn = 51;% 10 time steps required for accurate runaway solution - see rundke_dtn 0103 % graph.itdisp = 11:10:51;% display times - only used if opt.timevol > 0 0104 % 0105 % 0106 % 0107 opt.tnmin = 1e0;% > 0 for log scale 0108 % 0109 % opt.tn = [1e8,1e9];%[1e7,1e8];%1e10;% 10000 is time for asymptotic solution 0110 % opt.dtn = [161,90];%[141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0111 % graph.itdisp = [61:40:141,251];%[61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0 0112 % 0113 % opt.tn = [1e7,1e8];%1e10;% 10000 is time for asymptotic solution 0114 % opt.dtn = [141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0115 % graph.itdisp = [61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0 0116 % 0117 opt.tn = 1e4;%[1e7,1e8];%1e10;% 10000 is time for asymptotic solution 0118 opt.dtn = 81;%[141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0119 graph.itdisp = [21:20:81];%[61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0 0120 % 0121 % 0122 % 0123 % opt.tnmin = 1;%0;% > 0 for log scale 0124 % opt.tn = 1e8;%1e9;%1e10;% 10000 is time for asymptotic solution 0125 % opt.dtn = 17;%10*1i;%100*1i;%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0126 % graph.itdisp = [1,10];%[1,10,100];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0 0127 % 0128 % momentum space grid -> you can specify nmhu_S and either 0129 % - the momentum grid pn_S directly 0130 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0131 % 0132 grid.nmhu_S = 141;% number of pitch angle grid points 0133 grid.sfac1 = 80;%40;% % reduction factor for fine grid near abs(mhu) = 1 0134 grid.snmhu1 = 5; % number of fine grid points near abs(mhu) = 1 0135 % 0136 % momentum grid contruction option : 0137 % (0) do not adjust number of fine and coarse grid points 0138 % (1) keep fine grid step constant (typically, 0.2) 0139 % if isreal(np_tail) : keep coarse grid steps constant (typically, 1). Note : do not enforce np_S and np_tail variations 0140 % if isimag(np_tail), enforce continuity of grid step size at pnmax_S. Note : do not enforce np_S and np_tail variations 0141 % 0142 grid.opt = 1; 0143 % 0144 grid.np_S = 141;% number of bulk momentum grid points 0145 grid.pnmax_S = 28;%sqrt(((1 + Emax*1e3/mc2)^2 - 1)/betath2);% maximum of bulk momentum (normalized to pT) [pn = 28 corresponds to 1 MeV electrons for betath = 0.1] 0146 grid.np_tail = 4i;%960;%558;%331;%217;%103;%100;% number of tail momentum grid points 0147 grid.pnmax_S_tail = 1967;%988;%597;%401;205;% maximum of tail momentum (normalized to pT) [107,205,303,401,499,597,793,988,1184,1576,1967] corresponds to [5,10,15,20,25,30,40,50,60,80,100] MeV electrons for betath = 0.1 0148 % 0149 % avalanche modelling : 0150 % 0151 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator 0152 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated 0153 % - For pnmin0_KO = NaN, the limit is set to pnmax_S 0154 % 0155 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator 0156 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy' 0157 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term 0158 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO 0159 % 0160 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator 0161 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used 0162 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference) 0163 % 0164 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator 0165 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons. 0166 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO 0167 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO 0168 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number" 0169 % 0170 opt.avalanche_mode = false;%true;% 0171 opt.pnmin0_KO = NaN;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1 0172 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;% 0173 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator 0174 opt.pnmax2_KO = NaN;%100;% 0175 % 0176 path_simul = '';%Simulation path 0177 % 0178 if opt.bounce_mode == 0, 0179 id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect 0180 else 0181 id_equil = 'TScirc_e1';% toroidal equilibirum with a/R = 1 : then rho_S = r/R 0182 end 0183 id_wave = ''; 0184 % 0185 rho_S = 0.5; 0186 % 0187 opt.heavyions = true; 0188 opt.adjustB = true; 0189 % 0190 opt.keyboard = true; 0191 % 0192 graph.tnc = true; 0193 graph.xec = true; 0194 graph.xlog = 1; 0195 graph.xlim = [1e-3,1e2]; 0196 graph.xtick = 10.^(-3:2); 0197 % 0198 if alpha == 2, 0199 % 0200 % E/Ec = 2 0201 % 0202 graph.ylim_norm = [1e-25,1e1]; 0203 graph.ytick_norm = 10.^(-25:5:5); 0204 graph.ylim_gamma = [1e-26,1e-12]; 0205 graph.ytick_gamma = 10.^(-26:4:-12); 0206 graph.ylim_f = [1e-20,1e0];%[1e-45,1e5]; 0207 graph.ytick_f = 10.^(-20:4:0);%10.^(-45:5:5); 0208 elseif alpha == 2.5, 0209 % 0210 % E/Ec = 2.5 0211 % 0212 graph.ylim_norm = [1e-15,1e1]; 0213 graph.ytick_norm = 10.^(-15:3:0); 0214 graph.ylim_gamma = [1e-25,1e-10]; 0215 graph.ytick_gamma = 10.^(-25:3:-10); 0216 graph.ylim_f = [1e-45,1e5]; 0217 graph.ytick_f = 10.^(-45:5:5); 0218 elseif alpha == 3, 0219 % 0220 % E/Ec = 3 0221 % 0222 graph.ylim_norm = [1e-11,1e1];%[1e-25,1e1]; 0223 graph.ytick_norm = 10.^(-10:2:0);%10.^(-25:5:5); 0224 graph.ylim_gamma = [1e-12,1e-2]; 0225 graph.ytick_gamma = 10.^(-12:2:-2); 0226 graph.ylim_f = [1e-20,1e0];%[1e-35,1e5]; 0227 graph.ytick_f = 10.^(-20:4:0);%10.^(-35:5:5); 0228 % 0229 graph.dp_cyl = 0;%0.2; 0230 graph.cont = (0:3:100)*1e-11; 0231 % graph.cont = 1e-7*(-31:2:15)/8;%-3.875e-7:0.25e-7:1.875e-7;%-4e-7:0.25e-7:2e-7; 0232 graph.pnmax_disp = 300; 0233 graph.xlim_disp = [-5,20]; 0234 graph.ylim_disp = [0,10]; 0235 graph.xtick_disp = -5:5:20; 0236 graph.ytick_disp = 0:5:10; 0237 % 0238 elseif alpha == 5, 0239 % 0240 % E/Ec = 5 0241 % 0242 graph.ylim_norm = [1e-9,1e1]; 0243 graph.ytick_norm = 10.^(-10:2:2); 0244 graph.ylim_gamma = [1e-12,1e-4]; 0245 graph.ytick_gamma = 10.^(-12:2:-4); 0246 else 0247 graph.ylim_norm = NaN; 0248 graph.ylim_gamma = NaN; 0249 end 0250 % 0251 %************************************************************************************************************************************ 0252 % 0253 dp_bulk = grid.pnmax_S/(grid.np_S - 1); 0254 dp_fine = (grid.pnmax_S_tail - grid.pnmax_S)/grid.np_tail; 0255 pfac = imag(grid.np_tail); 0256 if grid.opt == 1, 0257 grid.np_S = 1 + round(grid.pnmax_S/dp_bulk); 0258 if pfac == 0, 0259 grid.np_tail = round((grid.pnmax_S_tail - grid.pnmax_S)/dp_fine); 0260 else 0261 grid.np_tail = round((grid.np_S - 1)*pfac*((grid.pnmax_S_tail/grid.pnmax_S)^(1/pfac) - 1)); 0262 grid.opt = pfac; 0263 end 0264 end 0265 % 0266 opt.pnmin_S = sqrt(((1 + Emin*1e3/mc2)^2 - 1)/betath2); 0267 wpr = -taurinv*1i;% real wpr : wpr, positive imaginary part : bhat, negative imaginary part : taurinv. 0268 % 0269 % RUN LUKE 0270 % 0271 [RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S); 0272 % 0273 %************************************************************************************************************************************ 0274 % 0275 % PROCESS DATA 0276 % 0277 fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul); 0278 % 0279 %************************************************************************************************************************************ 0280 % 0281 if p_opt > 0, 0282 filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat']; 0283 save(filename); 0284 info_dke_yp(2,['Data saved in ',filename]); 0285 end