rundke_runaway_taur_scan

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