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 % printing option 0022 % 0023 id_simul = 'Runaway_alpha3_sync';%Simulation ID 0024 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 0025 % 0026 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants 0027 % 0028 % calculation mode : (0) quasi steady state (1) time evolution 0029 % 0030 opt.timevol = 1; 0031 % 0032 % normalized electric field (with respect to the Critical field Ec) 0033 % 0034 alpha = 3;% E/Ec 0035 % 0036 % temperature betath2 = Te/mc2 0037 % (betath2 = 10^(-6) is validated for NR limit) 0038 % 0039 betath2 = 0.01;% Te/mc2 0040 % 0041 % the ratio of the syncrotron reaction force to the collisional drag scales as bhat^2 0042 % -> set bhat = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0043 % 0044 % note : 0045 % - without the syncrotron reaction force, and with a single cold ion 0046 % species of charge Zi, the normalized runaway rate - with respect to the 0047 % electron density and ee collision time - only depends upon Zi, alpha and betath2 0048 % 0049 % - with the syncrotron reaction force, an additional parameter is 0050 % provived by bhat 0051 % 0052 % - the link between wpr=wc2/wp2 and bhat is : 0053 % bhat^2 = wpr*(2/3)*betath^3/ln(Lambda) 0054 % 0055 % - the magnetic field being fixed, the density is adjusted in LUKE according to bhat 0056 % 0057 bhat = 0.03; 0058 % 0059 opt.synchro_mode = 1; 0060 % 0061 % ion charge (single species, cold ions) 0062 % 0063 Zi = 1; 0064 % 0065 Emin = NaN;%runaway energy threshold in MeV (limit where electrons are "counted") 0066 % 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] 0067 % 0068 % collision mode : 0069 % (0) : Relativistic Maxwellian background 0070 % (1) : High-velocity limit 0071 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0072 % 0073 opt.coll_mode = 0; 0074 opt.bounce_mode = 0;%1;% bounce-averaged calculation 0075 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0076 % 0077 % normalization mode : 0078 % (0) : "free" mode 0079 % (1) : "normalization of f0 maintained to its initial (maxwellian) value". DO NOT use in the runaway problem 0080 % (2) : "normalization of f0 maintained to its value at the previous time". Use only in the runaway problem with (+1i) 0081 % (... + 1i) : compensate for numerical error after each time step, acounting for boundary losses. 0082 opt.norm_mode_f = 0;%+1i; 0083 % 0084 % time grid (normalized to collision time) -> you can specify : 0085 % - linear arbitrary grid with tn array : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0 0086 % - linear arbitrary grid with dtn array : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10); 0087 % - linear grid with tnmax and grid step : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000; 0088 % - linear grid with tnmax and number of grid steps : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i; 0089 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10; 0090 % 0091 % opt.tnmin = 0;% > 0 for log scale 0092 % opt.tn = 10000;% 10000 is time for asymptotic solution 0093 % opt.dtn = 1000;% 10 time steps required for accurate runaway solution - see rundke_dtn 0094 % graph.itdisp = [1,5,10];% display times - only used if opt.timevol > 0 0095 % 0096 % opt.tnmin = 1;% > 0 for log scale 0097 % opt.tn = 1e5;% 10000 is time for asymptotic solution 0098 % opt.dtn = 51;% 10 time steps required for accurate runaway solution - see rundke_dtn 0099 % graph.itdisp = 11:10:51;% display times - only used if opt.timevol > 0 0100 % 0101 opt.tnmin = 1;% > 0 for log scale 0102 opt.tn = 1e6;% 10000 is time for asymptotic solution 0103 opt.dtn = 61;% 10 time steps required for accurate runaway solution - see rundke_dtn 0104 graph.itdisp = 11:10:61;% display times - only used if opt.timevol > 0 0105 % 0106 % momentum space grid -> you can specify nmhu_S and either 0107 % - the momentum grid pn_S directly 0108 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0109 % 0110 grid.nmhu_S = 101;% number of pitch angle grid points 0111 grid.np_S = 101;% number of bulk momentum grid points 0112 grid.pnmax_S = 20;%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] 0113 grid.np_tail = 100;%0;% number of tail momentum grid points 0114 grid.pnmax_S_tail = 20+100;% maximum of tail momentum (normalized to pT) [107,205,303,401,499,597] corresponds to [5,10,15,20,25,30] MeV electrons for betath = 0.1 0115 % 0116 % avalanche modelling : 0117 % 0118 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator 0119 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated 0120 % - For pnmin0_KO = NaN, the limit is set to pnmax_S 0121 % 0122 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator 0123 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy' 0124 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term 0125 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO 0126 % 0127 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator 0128 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used 0129 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference) 0130 % 0131 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator 0132 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons. 0133 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO 0134 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO 0135 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number" 0136 % 0137 opt.avalanche_mode = false;%true;% 0138 opt.pnmin0_KO = NaN;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1 0139 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;% 0140 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator 0141 opt.pnmax2_KO = NaN;%100;% 0142 % 0143 path_simul = '';%Simulation path 0144 % 0145 if opt.bounce_mode == 0, 0146 id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect 0147 else 0148 id_equil = 'TScirc_e1';% toroidal equilibirum with a/R = 1 : then rho_S = r/R 0149 end 0150 id_wave = ''; 0151 % 0152 rho_S = 0.5; 0153 % 0154 if alpha == 2, 0155 % 0156 % E/Ec = 2 0157 % 0158 graph.ylim_norm = [1e-25,1e1]; 0159 graph.ytick_norm = 10.^(-25:5:5); 0160 graph.ylim_gamma = [1e-26,1e-12]; 0161 graph.ytick_gamma = 10.^(-26:4:-12); 0162 graph.ylim_f = [1e-45,1e5]; 0163 graph.ytick_f = 10.^(-45:5:5); 0164 elseif alpha == 3, 0165 % 0166 % E/Ec = 3 0167 % 0168 graph.ylim_norm = [1e-25,1e1]; 0169 graph.ytick_norm = 10.^(-25:5:5); 0170 graph.ylim_gamma = [1e-32,1e-8]; 0171 graph.ytick_gamma = 10.^(-32:4:-8); 0172 graph.ylim_f = [1e-35,1e5]; 0173 graph.ytick_f = 10.^(-35:5:5); 0174 elseif alpha == 5, 0175 % 0176 % E/Ec = 5 0177 % 0178 graph.ylim_norm = [1e-9,1e1]; 0179 graph.ytick_norm = 10.^(-10:2:2); 0180 graph.ylim_gamma = [1e-12,1e-4]; 0181 graph.ytick_gamma = 10.^(-12:2:-4); 0182 else 0183 graph.ylim_norm = NaN; 0184 graph.ylim_gamma = NaN; 0185 end 0186 % 0187 %************************************************************************************************************************************ 0188 % 0189 opt.pnmin_S = sqrt(((1 + Emin*1e3/mc2)^2 - 1)/betath2); 0190 wpr = bhat*1i; 0191 % 0192 % RUN LUKE 0193 % 0194 [RR,mksa,Ec,f,tn,tRR,tnorm] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S); 0195 % 0196 %************************************************************************************************************************************ 0197 % 0198 % PROCESS DATA 0199 % 0200 fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul); 0201 % 0202 %************************************************************************************************************************************ 0203 % 0204 if p_opt > 0, 0205 filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat']; 0206 save(filename); 0207 info_dke_yp(2,['Data saved in ',filename]); 0208 end