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 permission = test_permissions_yp; 0013 % 0014 if ~permission 0015 disp('Please move the script to a local folder where you have write permission before to run it') 0016 return; 0017 end 0018 % 0019 % *********************** This part must be specified by the user ***************************** 0020 % 0021 scanstrs = {'alpha','sfac1','pnmax_S','np_tail','dtn','taurinv','betath2','zi'}; 0022 % 0023 % warning : 0024 % - do not vary np_tail with opt_grid = 1 0025 % - use taurinv, sfac1, pnmax_S, np_tail, dtn as scanstr1 only 0026 % - use betath2 as scanstr1 only if Emin is not NaN 0027 % 0028 scanstr1 = 'taurinv'; 0029 scanstr2 = 'alpha'; 0030 nan_enforce = [8,1;10,1];% artefacts in search for bump (linked to discretization) 0031 % 0032 % scanstr1 = 'zi'; 0033 % scanstr2 = 'alpha'; 0034 % nan_enforce = [11,2];%NaN(0,2);% artefacts in search for bumb (linked to discretization) 0035 % 0036 % MDCE 0037 % 0038 mdce_mode = -1;%0;%3,% 0039 remnum = 0; 0040 returnmode = 0; 0041 transpmode = true; 0042 nproc_enforce = 4; 0043 % 0044 % scan parameters 0045 % 0046 scanvar.dtn = 'opt.dtn'; 0047 scanval.dtn = 1 + 8*[1:4,8:4:40]; 0048 scanlab.dtn = '\Deltat/\tau_c'; 0049 xlim.dtn = [0,350]; 0050 xlog.dtn = 0; 0051 xtick.dtn = 0:50:350; 0052 % 0053 scanvar.np_tail = 'grid.np_tail'; 0054 scanval.np_tail = (988-28)./2.^(0:6); 0055 scanlab.np_tail = 'n_p (tail)'; 0056 xlim.np_tail = (988-28)./2.^[7,-1]; 0057 xlog.np_tail = 1; 0058 xtick.np_tail = fliplr(scanval.np_tail); 0059 % 0060 scanvar.pnmax_S = 'grid.pnmax_S'; 0061 scanval.pnmax_S = 7:7:70;%140; 0062 scanlab.pnmax_S = 'p_{max}/p_T (fine grid)'; 0063 xlim .pnmax_S= [0,77];%[0,147];% 0064 xlog.pnmax_S = 0; 0065 xtick.pnmax_S = scanval.pnmax_S(2:2:end); 0066 % 0067 scanvar.sfac1 = 'grid.sfac1'; 0068 scanval.sfac1 = [1,2,4,6,10,20,40,60,100]; 0069 scanlab.sfac1 = '\Delta\xi_{fine}/\Delta\xi_0'; 0070 xlim.sfac1 = [0.5,200];%[0,147];% 0071 xlog.sfac1 = 1; 0072 xtick.sfac1 = scanval.sfac1; 0073 % 0074 scanvar.alpha = 'alpha'; 0075 scanval.alpha = (1.5:0.25:4);%(1.5:0.5:6);% 0076 scanlab.alpha = 'E/E_c'; 0077 scanind.alpha = 2; 0078 xlim.alpha = [1.25,4.25];%[1,7];% 0079 xlog.alpha = 0; 0080 xtick.alpha = scanval.alpha(1:2:11); 0081 iscandisp.alpha = 1:2:11; 0082 % 0083 scanvar.betath2 = 'betath2'; 0084 scanval.betath2 = 10.^(-3:0.25:-1); 0085 scanlab.betath2 = 'E/E_c'; 0086 scanind.betath2 = 3; 0087 xlim.betath2 = 10.^[-3.25,-0.75]; 0088 xlog.betath2 = 1; 0089 xtick.betath2 = scanval.betath2(1:2:9); 0090 % 0091 scanvar.taurinv = 'taurinv'; 0092 scanval.taurinv = 0.3:0.15:1.95;%0.3:0.1:1.5; 0093 scanlab.taurinv = '\tau_r^{-1}'; 0094 xlim.taurinv = [0.25,2.0];%[0.2,1.6]; 0095 xlog.taurinv = 0; 0096 xtick.taurinv = scanval.taurinv(1:2:11);%scanval.taurinv(1:3:13);% 0097 iscandisp.taurinv = 1:2:11; 0098 % 0099 scanvar.zi = 'Zi'; 0100 scanval.zi = 1:0.2:3;%1:0.25:4;%1:10;% 0101 scanlab.zi = 'Z_i'; 0102 scanind.zi = 5; 0103 xlim.zi = [0.9,3.1];%[0.5,4.5];%[0,11];% 0104 xlog.zi = 0; 0105 xtick.zi =1:0.4:3;%0.5:0.5:4.5; 0106 iscandisp.zi = 1:2:11; 0107 % 0108 graph.cont_norm = 10.^(-24:0); 0109 graph.cont_norm_lab = 10.^(-24:4:0); 0110 graph.cont_gamma = 10.^(-30:2:0); 0111 graph.cont_gamma_lab = 10.^(-30:5:0); 0112 % 0113 % printing option 0114 % 0115 id_simul0 = 'Runaway_scan';%Simulation ID 0116 p_opt = -1;%2;% printing and saving options : (-1) do nothing (0) print figures (1) print figures and save figures and results (2) save figures and results 0117 % 0118 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants 0119 % 0120 % calculation mode : (0) quasi steady state (1) time evolution 0121 % 0122 opt.timevol = 1; 0123 % 0124 % normalized electric field (with respect to the Critical field Ec) 0125 % 0126 alpha = 3;% E/Ec 0127 % 0128 % temperature betath2 = Te/mc2 0129 % (betath2 = 10^(-6) is validated for NR limit) 0130 % 0131 betath2 = 0.01;% Te/mc2 0132 % 0133 % the ratio of the syncrotron reaction force to the collisional drag scales 0134 % as taurinv = 1/taur 0135 % -> set taurinv = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0136 % 0137 % note : 0138 % - without the syncrotron reaction force, and with a single cold ion 0139 % species of charge Zi, the normalized runaway rate - with respect to the 0140 % electron density and ee collision time - only depends upon Zi, alpha and betath2 0141 % 0142 % - with the syncrotron reaction force, an additional parameter is 0143 % provived by taurinv 0144 % 0145 % - the link between wpr=wc2/wp2 and taurinv is : 0146 % taurinv = wpr*(2/3)/ln(Lambda) 0147 % 0148 % - the magnetic field being fixed, the density is adjusted in LUKE according to bhat 0149 % 0150 taurinv = 0.6;%0.1;% 0151 % 0152 opt.synchro_mode = 1; 0153 % 0154 % ion charge (single species, cold ions) 0155 % 0156 Zi = 1; 0157 % 0158 Emin = NaN;%runaway energy threshold in MeV (limit where electrons are "counted") 0159 % 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] 0160 % 0161 % collision mode : 0162 % (0) : Relativistic Maxwellian background 0163 % (1) : High-velocity limit 0164 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0165 % 0166 opt.coll_mode = 0; 0167 opt.bounce_mode = 0;%1;% bounce-averaged calculation 0168 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0169 % 0170 % normalization mode : 0171 % (0) : "free" mode 0172 % (1) : "normalization of f0 maintained to its initial (maxwellian) value". DO NOT use in the runaway problem 0173 % (2) : "normalization of f0 maintained to its value at the previous time". Use only in the runaway problem with (+1i) 0174 % (... + 1i) : compensate for numerical error after each time step, acounting for boundary losses. 0175 opt.norm_mode_f = +1i; 0176 % 0177 % time grid (normalized to collision time) -> you can specify : 0178 % - linear arbitrary grid with tn array : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0 0179 % - linear arbitrary grid with dtn array : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10); 0180 % - linear grid with tnmax and grid step : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000; 0181 % - linear grid with tnmax and number of grid steps : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i; 0182 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10; 0183 % 0184 % opt.tnmin = 0;% > 0 for log scale 0185 % opt.tn = 10000;% 10000 is time for asymptotic solution 0186 % opt.dtn = 1000;% 10 time steps required for accurate runaway solution - see rundke_dtn 0187 % 0188 % opt.tnmin = 1;% > 0 for log scale 0189 % opt.tn = 1e5;% 10000 is time for asymptotic solution 0190 % opt.dtn = 51;% 10 time steps required for accurate runaway solution - see rundke_dtn 0191 % 0192 opt.tnmin = 1;% > 0 for log scale 0193 opt.tn = [1e7,1e8];%1e10;% 10000 is time for asymptotic solution 0194 opt.dtn = [141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0195 % 0196 % opt.tnmin = 0;% > 0 for log scale 0197 % opt.tn = 1e8;%1e10;% 10000 is time for asymptotic solution 0198 % opt.dtn = 10*1i;%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn 0199 % 0200 % momentum space grid -> you can specify nmhu_S and either 0201 % - the momentum grid pn_S directly 0202 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0203 % 0204 grid.nmhu_S = 141;% number of pitch angle grid points 0205 grid.sfac1 = 80;%40; % reduction factor for fine grid near abs(mhu) = 1 0206 grid.snmhu1 = 5; % number of fine grid points near abs(mhu) = 1 0207 % 0208 % momentum grid contruction option : 0209 % (0) do not adjust number of fine and coarse grid points 0210 % (1) keep fine grid step constant (typically, 0.2) 0211 % if isreal(np_tail) : keep coarse grid steps constant (typically, 1). Note : do not enforce np_S and np_tail variations 0212 % if isimag(np_tail), enforce continuity of grid step size at pnmax_S. Note : do not enforce np_S and np_tail variations 0213 % 0214 grid.opt = 1; 0215 % 0216 grid.np_S = 141;% number of bulk momentum grid points 0217 grid.pnmax_S = 28;%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] 0218 grid.np_tail = 4i;%960;%240;%558;%331;%217;%103;%100;% number of tail momentum grid points 0219 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 0220 % 0221 % avalanche modelling : 0222 % 0223 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator 0224 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated 0225 % - For pnmin0_KO = NaN, the limit is set to pnmax_S 0226 % 0227 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator 0228 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy' 0229 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term 0230 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO 0231 % 0232 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator 0233 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used 0234 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference) 0235 % 0236 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator 0237 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons. 0238 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO 0239 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO 0240 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number" 0241 % 0242 opt.avalanche_mode = false;%true;% 0243 opt.pnmin0_KO = NaN;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1 0244 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;% 0245 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator 0246 opt.pnmax2_KO = NaN;%100;% 0247 % 0248 path_simul = '';%Simulation path 0249 % 0250 if opt.bounce_mode == 0, 0251 id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect 0252 else 0253 id_equil = 'TScirc_e1';% toroidal equilibirum with a/R = 1 : then rho_S = r/R 0254 end 0255 id_wave = ''; 0256 % 0257 rho_S = 0.5; 0258 % 0259 opt.heavyions = true; 0260 opt.adjustB = true; 0261 % 0262 %************************************************************************************************************************************ 0263 % 0264 nscan1 = length(scanval.(scanstr1)); 0265 nscan2 = length(scanval.(scanstr2)); 0266 id_simul_scan = [id_simul0,'_',scanstr1,'_',scanstr2]; 0267 % 0268 RR = cell(nscan1,nscan2); 0269 mksa = cell(nscan1,nscan2); 0270 Ec = cell(nscan1,nscan2); 0271 f = cell(nscan1,nscan2); 0272 tn = cell(nscan1,nscan2); 0273 tRR = cell(nscan1,nscan2); 0274 tnorm = cell(nscan1,nscan2); 0275 fgrid = cell(nscan1,nscan2); 0276 % 0277 dp_bulk = grid.pnmax_S/(grid.np_S - 1); 0278 dp_fine = (grid.pnmax_S_tail - grid.pnmax_S)/grid.np_tail; 0279 pfac = imag(grid.np_tail); 0280 % 0281 for iscan1 = 1:nscan1, 0282 evalstr = [scanvar.(scanstr1),' = ',num2str(scanval.(scanstr1)(iscan1)),';']; 0283 disp(['Scan ',num2str(iscan1),'/',num2str(nscan1),' : ',evalstr]) 0284 eval(evalstr); 0285 % 0286 if grid.opt == 1, 0287 grid.np_S = 1 + round(grid.pnmax_S/dp_bulk); 0288 if pfac == 0, 0289 grid.np_tail = round((grid.pnmax_S_tail - grid.pnmax_S)/dp_fine); 0290 else 0291 grid.np_tail = round((grid.np_S - 1)*pfac*((grid.pnmax_S_tail/grid.pnmax_S)^(1/pfac) - 1)); 0292 grid.opt = pfac; 0293 end 0294 end 0295 % 0296 id_simul = [id_simul_scan,'_',num2str(scanval.(scanstr1)(iscan1))]; 0297 % 0298 opt.pnmin_S = sqrt(((1 + Emin*1e3/mc2)^2 - 1)/betath2); 0299 wpr = -taurinv*1i;% real wpr : wpr, positive imaginary part : bhat, negative imaginary part : taurinv. 0300 % 0301 % RUN LUKE 0302 % 0303 dkepath = lukeschedulerparam('',mdce_mode,remnum,returnmode,transpmode); 0304 % 0305 if ~isnan(nproc_enforce), 0306 for irem = 1:length(dkepath.remote), 0307 dkepath.remote(irem).nproc = nproc_enforce; 0308 end 0309 end 0310 % 0311 dkecluster = clustermode_luke('','run_lukert',dkepath); 0312 % 0313 [flag,RR(iscan1,:),mksa(iscan1,:),Ec(iscan1,:),f(iscan1,:),tn(iscan1,:),tRR(iscan1,:),tnorm(iscan1,:),fgrid(iscan1,:)] = ... 0314 mdce_run('frundke_runaway',{id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S},scanind.(scanstr2),scanval.(scanstr2),dkecluster); 0315 % 0316 end 0317 % 0318 %************************************************************************************************************************************ 0319 % 0320 % PROCESS DATA 0321 % 0322 fproc_runaway_scan_2D(scanstr1,scanstr2,scanvar,scanval,scanlab,RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,iscandisp,xlim,xlog,xtick,p_opt,id_simul_scan,nan_enforce); 0323 % 0324 %************************************************************************************************************************************ 0325 % 0326 if p_opt > 0, 0327 % 0328 for iscan1 = 1:nscan1, 0329 for iscan2 = 1:nscan2, 0330 f{iscan1,iscan2} = f{iscan1,iscan2}(end);% to save space 0331 end 0332 end 0333 % 0334 filename = [path_simul,'DKE_RESULTS_',id_simul_scan,'.mat']; 0335 save(filename); 0336 info_dke_yp(2,['Data saved in ',filename]); 0337 end