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 % normalized electric field (with respect to the Critical field Ec) 0027 % 0028 alpha = [10]*1e4;%E/Ec 0029 % 0030 % temperature betath2 = Te/mc2 0031 % (betath2 = 10^(-6) is validated for NR limit) 0032 % 0033 betath2 = 1e-6; 0034 % 0035 % ion charge (single species, cold ions) 0036 % 0037 Zi = 2;%3;% 0038 % 0039 % collision mode : 0040 % (0) : Relativistic Maxwellian background 0041 % (1) : High-velocity limit 0042 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0043 % 0044 coll_mode = 0; 0045 % 0046 p_opt = -1; 0047 % 0048 % ***********************This part must be specified by the user, run make files if necessary) ***************************** 0049 % 0050 id_simul = 'Runaway_critical';%Simulation ID 0051 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0052 % 0053 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty) 0054 rho_S = [0.5];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty) 0055 % 0056 id_path = '';%For all paths used by DKE solver 0057 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0058 % 0059 id_equil = 'TScyl';%For plasma equilibrium 0060 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0061 % 0062 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters 0063 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0064 % 0065 id_display = 'PARTIAL_VISUAL';%For output code display 0066 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0067 % 0068 id_ohm = '';%For Ohmic electric contribution 0069 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0070 % 0071 ids_wave = {''};%For RF waves contribution (put all the type of waves needed) 0072 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0073 % 0074 id_transpfaste = '';%For fast electron radial transport 0075 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0076 % 0077 id_ripple = '';%For fast electron magnetic ripple losses 0078 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0079 % 0080 %************************************************************************************************************************************ 0081 %************************************************************************************************************************************ 0082 %************************************************************************************************************************************ 0083 % 0084 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple); 0085 % 0086 %************************************************************************** 0087 %************************************************************************** 0088 % 0089 if exist('dmumpsmex');dkeparam.invproc = -2;end 0090 % 0091 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0092 dkeparam.norm_mode_f = 1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes 0093 dkeparam.tn = 10000;%time for asymptotic solution with norm_mode_f = 1 0094 dkeparam.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn 0095 dkeparam.coll_mode = coll_mode;% mode for the collision operator 0096 % 0097 dkeparam.psin_S = psin_S; 0098 dkeparam.rho_S = rho_S; 0099 % 0100 equil.pzTi = 1e-10*ones(size(equil.pzTi));% cold ions 0101 % 0102 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants 0103 % 0104 equil.pTe = betath2*mc2*ones(size(equil.pTe)); 0105 % 0106 equil.zZi = [1,1,1,Zi]; 0107 equil.pzni = [0;0;0;1/Zi]*equil.pne;% single ion species 0108 % 0109 epsi = alpha*betath2; 0110 ohm = ohm_flat(equil,epsi); 0111 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0112 % 0113 RR = dke_out.XxRR_fsav(end,:); 0114 % 0115 epsi_kulsrud = [0.04;0.06;0.08;0.10]; 0116 Zi_kulsrud =[1,2,3,10]; 0117 RR_kulsrud = 1e-4*[ 0.01914,NaN ,NaN ,NaN ;...%Kulsrud (PRL, 31,11, (1972) 690) 0118 0.5411 ,0.2611,NaN ,NaN ;... 0119 3.177 ,1.735 ,1.047,0.09 ;... 0120 10.04 ,5.839 ,3.757,0.449]; 0121 % 0122 iepsi = find(abs((epsi - epsi_kulsrud)./epsi_kulsrud) < 10*eps); 0123 % 0124 if coll_mode == 0 && abs((betath2 - 1e-6)/1e-6) < 10*eps && ~isempty(iepsi) && any(Zi == Zi_kulsrud), 0125 RR_k = RR_kulsrud(iepsi,Zi == Zi_kulsrud); 0126 else 0127 RR_k = ''; 0128 end 0129 % 0130 % additional parameters 0131 % 0132 % Dreicer field ED = Ec/betath2 0133 % Critical momentum at mhu=1 : 0134 % - NR limit : pc = 1/sqrt(betath2*alpha) 0135 % - Relativistic : pc = 1/sqrt(betath2*alpha - betath2)) 0136 % 0137 %************************************************************************************************************************************ 0138 % 0139 format 0140 % 0141 delete res_critical 0142 % 0143 diary res_critical 0144 % 0145 disp('Simulation parameters : ') 0146 disp(['coll_mode = ',num2str(coll_mode),' ; Te/mc2 = ',num2str(betath2),' ; E/Ec = ',num2str(alpha),' ; E/ED = ',num2str(epsi),' ; Zi = ',num2str(Zi)]) 0147 disp(['NR limit : pc/pT = ',num2str(1/sqrt(epsi)),' ; pmax/pc = ',num2str(dkeparam.pnmax_S*sqrt(epsi))])%NR 0148 disp(['Relativ. : pc/pT = ',num2str(1/sqrt(epsi - betath2)),' ; pmax/pc = ',num2str(dkeparam.pnmax_S*sqrt(epsi - betath2))])%rel 0149 disp(' ') 0150 disp(['--> LUKE RR : ',num2str(RR)]) 0151 if ~isempty(RR_k), 0152 disp(['--> Kulsrud RR : ',num2str(RR_k)]) 0153 end 0154 disp('----------------------') 0155 % 0156 diary off 0157 % 0158 %************************************************************************************************************************************ 0159 % 0160 filename = [path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']; 0161 save(filename); 0162 info_dke_yp(2,['Data saved in ',filename]);