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_kulsrud_time';%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 = 6e4;%1e5;%E/Ec 0035 % 0036 % temperature betath2 = Te/mc2 0037 % (betath2 = 10^(-6) is validated for NR limit) 0038 % 0039 betath2 = 1e-6; 0040 % 0041 % the ratio of the syncrotron reaction force to the collisional drag scales as wc2/wp2 0042 % -> set wpr = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0043 % 0044 wpr = 0;%([Bt]*qe/me)^2/([ne19]*1e19*qe^2/(me*e0));%wc2/wp2 0045 % 0046 opt.synchro_mode = 0; 0047 % 0048 % ion charge (single species, cold ions) 0049 % 0050 Zi = 1;%2;% 0051 % 0052 % collision mode : 0053 % (0) : Relativistic Maxwellian background 0054 % (1) : High-velocity limit 0055 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0056 % 0057 opt.coll_mode = 0; 0058 opt.bounce_mode = 0;%1;% bounce-averaged calculation 0059 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0060 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 0061 % 0062 % time grid (normalized to collision time) -> you can specify : 0063 % - linear arbitrary grid with tn array : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0 0064 % - linear arbitrary grid with dtn array : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10); 0065 % - linear grid with tnmax and grid step : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000; 0066 % - linear grid with tnmax and number of grid steps : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i; 0067 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10; 0068 % 0069 opt.tnmin = 1;%0;% > 0 for log scale 0070 opt.tn = 1e4;%10000;% 10000 is time for asymptotic solution 0071 opt.dtn = 81;%1000;% 10 time steps required for accurate runaway solution - see rundke_dtn 0072 % 0073 % momentum space grid -> you can specify nmhu_S and either 0074 % - the momentum grid pn_S directly 0075 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0076 % 0077 grid.nmhu_S = 101;% number of pitch angle grid points 0078 grid.np_S = 101;% number of bulk momentum grid points 0079 grid.pnmax_S = 20;% maximum of bulk momentum (normalized to pT) 0080 grid.np_tail = 0;%100;% number of tail momentum grid points 0081 grid.pnmax_S_tail = 5i;% maximum of tail momentum (normalized to pT) 0082 % 0083 path_simul = '';%Simulation path 0084 % 0085 id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect 0086 id_wave = ''; 0087 % 0088 rho_S = 0.5; 0089 % 0090 % graph.itdisp = 20:20:60;%[1,5,10];% display times - only used if opt.timevol > 0 0091 graph.itdisp = 21:20:81;%[1,5,10];% display times - only used if opt.timevol > 0 0092 % 0093 if alpha == 1e5 && Zi == 2, 0094 graph.ylim_norm = [1e-4,2e-0]; 0095 graph.ytick_norm = 10.^(-4:1); 0096 % 0097 graph.ylim_gamma = [5e-5,2e-3]; 0098 graph.ytick_gamma = 10.^(-5:-2); 0099 elseif alpha == 1e5 && Zi == 2, 0100 graph.ylim_norm = [1e-4,2e-0]; 0101 graph.ytick_norm = 10.^(-4:1); 0102 % 0103 graph.ylim_gamma = [1e-6,1e-4]; 0104 graph.ytick_gamma = 10.^(-6:-4); 0105 else 0106 graph.ylim_norm = NaN; 0107 graph.ylim_gamma = NaN; 0108 end 0109 % 0110 %************************************************************************************************************************************ 0111 % 0112 % RUN LUKE 0113 % 0114 [RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,ZXXS] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S); 0115 % 0116 %************************************************************************************************************************************ 0117 % 0118 % PROCESS DATA 0119 % 0120 fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,ZXXS,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul); 0121 % 0122 %************************************************************************************************************************************ 0123 % 0124 if p_opt > 0, 0125 filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat']; 0126 save(filename); 0127 info_dke_yp(2,['Data saved in ',filename]); 0128 end