This script is used to run LUKE for the runaway problem by J. Decker, Y. Peysson and E. Nilsson Note : A comparison with Kulsrud is provided for the follwoing parameters : - betath2 = 1e-6 (NR limit) - alpha among [4,6,8,10]*1e4; - Zi among [1,2,3,10] - coll_mode = 0 (Relativistic Maxwellian background)
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 % Note : A comparison with Kulsrud is provided for the follwoing parameters : 0007 % 0008 % - betath2 = 1e-6 (NR limit) 0009 % - alpha among [4,6,8,10]*1e4; 0010 % - Zi among [1,2,3,10] 0011 % - coll_mode = 0 (Relativistic Maxwellian background) 0012 % 0013 clear all 0014 clear mex 0015 clear functions 0016 close all 0017 warning('off') 0018 % 0019 permission = test_permissions_yp; 0020 % 0021 if ~permission 0022 disp('Please move the script to a local folder where you have write permission before to run it') 0023 return; 0024 end 0025 % 0026 % *********************** This part must be specified by the user ***************************** 0027 % 0028 % 0029 % normalized electric field (with respect to the Critical field Ec) 0030 % 0031 alpha = 1e5;%E/Ec 0032 % 0033 % temperature betath2 = Te/mc2 0034 % (betath2 = 10^(-6) is validated for NR limit) 0035 % 0036 betath2 = 1e-6;% 0037 % 0038 % the ratio of the syncrotron reaction force to the collisional drag scales as wc2/wp2 0039 % -> set wpr = 0 to ignore syncrotron reaction force 0040 % 0041 wpr = 0;%wc2/wp2 0042 % 0043 % ion charge (single species, cold ions) 0044 % 0045 Zi = 2;% 0046 % 0047 % collision mode : 0048 % (0) : Relativistic Maxwellian background 0049 % (1) : High-velocity limit 0050 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0051 % 0052 opt.coll_mode = 0; 0053 opt.bounce_mode = 0;% bounce-averaged calculation 0054 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0055 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 0056 % 0057 opt.tn = 10000;%10000 is time for asymptotic solution 0058 opt.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn 0059 % 0060 id_simul = 'Runaway_kulsrud';%Simulation ID 0061 path_simul = '';%Simulation path 0062 % 0063 %************************************************************************************************************************************ 0064 % 0065 RR = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt); 0066 % 0067 %************************************************************************************************************************************ 0068 % 0069 epsi_kulsrud = [0.04;0.06;0.08;0.10]; 0070 Zi_kulsrud =[1,2,3,10]; 0071 RR_kulsrud = 1e-4*[ 0.01914,NaN ,NaN ,NaN ;...%Kulsrud (PRL, 31,11, (1972) 690) 0072 0.5411 ,0.2611,NaN ,NaN ;... 0073 3.177 ,1.735 ,1.047,0.09 ;... 0074 10.04 ,5.839 ,3.757,0.449]; 0075 % 0076 iepsi = find(abs((alpha*betath2 - epsi_kulsrud)./epsi_kulsrud) < 10*eps); 0077 % 0078 if opt.coll_mode == 0 && opt.bounce_mode == 0 && abs((betath2 - 1e-6)/1e-6) < 10*eps && ~isempty(iepsi) && any(Zi == Zi_kulsrud), 0079 RR_k = RR_kulsrud(iepsi,Zi == Zi_kulsrud); 0080 else 0081 RR_k = ''; 0082 end 0083 % 0084 % additional info 0085 % 0086 % Dreicer field ED = Ec/betath2 0087 % Critical momentum at mhu=1 : 0088 % - NR limit : pc = 1/sqrt(betath2*alpha) 0089 % - Relativistic : pc = 1/sqrt(betath2*alpha - betath2)) 0090 % 0091 %************************************************************************************************************************************ 0092 % 0093 format 0094 % 0095 delete res_runaway_kulsrud 0096 % 0097 diary res_runaway_kulsrud 0098 % 0099 disp('Simulation parameters : ') 0100 disp(['coll_mode = ',num2str(opt.coll_mode),' ; bounce_mode = ',num2str(opt.bounce_mode)]) 0101 disp(['Te/mc2 = ',num2str(betath2),' ; E/Ec = ',num2str(alpha),' ; E/ED = ',num2str(alpha*betath2),' ; Zi = ',num2str(Zi)]) 0102 disp(['NR limit : pc/pT = ',num2str(1/sqrt(alpha*betath2))])%NR 0103 disp(['Relativ. : pc/pT = ',num2str(1/sqrt(alpha*betath2 - betath2))])%rel 0104 disp(' ') 0105 disp(['--> LUKE RR : ',num2str(RR)]) 0106 if ~isempty(RR_k), 0107 disp(['--> Kulsrud RR : ',num2str(RR_k)]) 0108 end 0109 disp('----------------------') 0110 % 0111 diary off 0112 % 0113 %************************************************************************************************************************************ 0114 %