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_test_a5_lin_fine_long';%Simulation ID 0024 p_opt = 2;% 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 Emin = 1;%runaway energy threshold in MeV 0029 Emax = 1;%grid energy threshold in MeV 0030 % 0031 % calculation mode : (0) quasi steady state (1) time evolution 0032 % 0033 opt.timevol = 1; 0034 % 0035 % normalized electric field (with respect to the Critical field Ec) 0036 % 0037 alpha = 5;%E/Ec 0038 % 0039 % temperature betath2 = Te/mc2 0040 % (betath2 = 10^(-6) is validated for NR limit) 0041 % 0042 betath2 = 0.01; 0043 % 0044 % the ratio of the syncrotron reaction force to the collisional drag scales as wc2/wp2 0045 % -> set wpr = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0046 % 0047 wpr = 0;%([Bt]*qe/me)^2/([ne19]*1e19*qe^2/(me*e0));%wc2/wp2 0048 % 0049 opt.synchro_mode = 0; 0050 % 0051 % ion charge (single species, cold ions) 0052 % 0053 Zi = 1; 0054 % 0055 % collision mode : 0056 % (0) : Relativistic Maxwellian background 0057 % (1) : High-velocity limit 0058 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0059 % 0060 opt.coll_mode = 0; 0061 opt.bounce_mode = 0;% bounce-averaged calculation 0062 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0063 opt.norm_mode_f = 0;%Local normalization of f0 at each iteration (0) no, the default value when the numerical conservative scheme is correct, (1) yes 0064 % 0065 % time grid (normalized to collision time) -> you can specify : 0066 % - linear arbitrary grid with tn array : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0 0067 % - linear arbitrary grid with dtn array : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10); 0068 % - linear grid with tnmax and grid step : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000; 0069 % - linear grid with tnmax and number of grid steps : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i; 0070 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10; 0071 % 0072 % opt.tnmin = 1;%0;% > 0 for log scale 0073 % opt.tn = 1e4;%1e5;% 10000 is time for asymptotic solution 0074 % opt.dtn = 81;%1001;% 10 time steps required for accurate runaway solution - see rundke_dtn 0075 % graph.itdisp = 21:20:81;%201:200:1001;% display times - only used if opt.timevol > 0 0076 % 0077 opt.tnmin = 0;% > 0 for log scale 0078 opt.tn = 10000;%100000;% 10000 is time for asymptotic solution 0079 opt.dtn = 10;%1000;%1000;% 10 time steps required for accurate runaway solution - see rundke_dtn 0080 graph.itdisp = 200:200:1000;%20:20:100;%[1,5,10];% display times - only used if opt.timevol > 0 0081 % 0082 % momentum space grid -> you can specify nmhu_S and either 0083 % - the momentum grid pn_S directly 0084 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0085 % 0086 % grid.nmhu_S = 101;% number of pitch angle grid points 0087 % grid.np_S = 101;% number of bulk momentum grid points 0088 % %grid.pnmax_S = 28;% maximum of bulk momentum (normalized to pT) [pn = 28 corresponds to 1 MeV electrons for betath = 0.1] 0089 % grid.pnmax_S = sqrt(((1 + Emax*1e3/mc2)^2 - 1)/betath2);% maximum of bulk momentum (normalized to pT) 0090 % grid.np_tail = 0;%100;% number of tail momentum grid points 0091 % grid.pnmax_S_tail = 5i;% maximum of tail momentum (normalized to pT) 0092 % 0093 grid.nmhu_S = 101;% number of pitch angle grid points 0094 grid.np_S = 141;% number of bulk momentum grid points 0095 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] 0096 grid.np_tail = 331;%100;% number of tail momentum grid points 0097 grid.pnmax_S_tail = 597;% 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 0098 % 0099 opt.pnmin_S = grid.pnmax_S;% to only count RE beyond 1MeV 0100 % 0101 % avalanche modelling : 0102 % 0103 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator 0104 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated 0105 % - For pnmin0_KO = NaN, the limit is set to pnmax_S 0106 % 0107 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator 0108 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy' 0109 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term 0110 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO 0111 % 0112 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator 0113 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used 0114 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference) 0115 % 0116 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator 0117 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons. 0118 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO 0119 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO 0120 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number" 0121 % 0122 opt.avalanche_mode = false;%true;% 0123 opt.pnmin0_KO = 28;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1 0124 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;% 0125 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator 0126 opt.pnmax2_KO = NaN;%100;% 0127 % 0128 path_simul = '';%Simulation path 0129 % 0130 id_equil = '';%'TScirc';%'TScyl';% 0131 id_wave = ''; 0132 % 0133 rho_S = 0.005; 0134 % 0135 if alpha == 10, 0136 % 0137 % E/Ec = 10 0138 % 0139 graph.ylim_norm = [1e-4,2e-0]; 0140 graph.ytick_norm = 10.^(-4:1); 0141 % 0142 graph.ylim_gamma = [1e-4,5e-3]; 0143 graph.ytick_gamma = 10.^(-5:-2); 0144 % 0145 elseif alpha == 5, 0146 % 0147 % E/Ec = 5 0148 % 0149 graph.ylim_norm = [1e-9,1e1]; 0150 graph.ytick_norm = 10.^(-10:2:2); 0151 % 0152 graph.ylim_gamma = [1e-12,1e-4]; 0153 graph.ytick_gamma = 10.^(-12:2:-4); 0154 else 0155 graph.ylim_norm = NaN; 0156 graph.ylim_gamma = NaN; 0157 end 0158 % 0159 %************************************************************************************************************************************ 0160 % 0161 % RUN LUKE 0162 % 0163 [RR,mksa,Ec,f,tn,tRR,tnorm] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S); 0164 % 0165 %************************************************************************************************************************************ 0166 % 0167 % PROCESS DATA 0168 % 0169 %fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul); 0170 fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,'','','','',alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul) 0171 % 0172 %************************************************************************************************************************************ 0173 % 0174 if p_opt > 0, 0175 filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat']; 0176 save(filename); 0177 info_dke_yp(2,['Data saved in ',filename]); 0178 end