rundke_runaway_taur_streamtest

PURPOSE ^

This script is used to run LUKE for the runaway problem

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 This script is used to run LUKE for the runaway problem

 by J. Decker, Y. Peysson and E. Nilsson

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 % *********************** This part must be specified by the user *****************************
0013 %
0014 % printing option
0015 %
0016 permission = test_permissions_yp;
0017 %
0018 if ~permission 
0019     disp('Please move the script to a local folder where you have write permission before to run it')
0020     return;
0021 end
0022 %
0023 id_simul = 'Runaway_alpha3_nosync';%_coll2%Simulation ID
0024 % id_simul = 'Runaway_alpha3_sync_0p6';%_coll2_norm1_Zi4%Simulation ID
0025 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
0026 %
0027 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants
0028 %
0029 % calculation mode : (0) quasi steady state (1) time evolution
0030 %
0031 opt.timevol = 1;
0032 %
0033 % normalized electric field (with respect to the Critical field Ec)
0034 %
0035 alpha = 3;%3.25;% E/Ec
0036 %
0037 % temperature betath2 = Te/mc2
0038 % (betath2 = 10^(-6) is validated for NR limit)
0039 %
0040 betath2 = 0.01;% Te/mc2
0041 %
0042 % the ratio of the syncrotron reaction force to the collisional drag scales
0043 % as taurinv = 1/taur
0044 % -> set taurinv = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force
0045 %
0046 % note :
0047 %   - without the syncrotron reaction force, and with a single cold ion
0048 %   species of charge Zi, the normalized runaway rate - with respect to the
0049 %   electron density and ee collision time - only depends upon Zi, alpha and betath2
0050 %
0051 %   - with the syncrotron reaction force, an additional parameter is
0052 %   provived by taurinv
0053 %
0054 %   - the link between wpr=wc2/wp2 and taurinv is :
0055 %       taurinv = wpr*(2/3)/ln(Lambda)
0056 %
0057 %   - the magnetic field being fixed, the density is adjusted in LUKE according to bhat
0058 %
0059 taurinv = 0.6;%0.3;%1;%0.7523;%0.9;%
0060 %
0061 opt.synchro_mode = 0;
0062 % opt.synchro_mode = 1;
0063 %
0064 % ion charge (single species, cold ions)
0065 %
0066 Zi = 1;%4;%
0067 %
0068 Emin = NaN;%runaway energy threshold in MeV (limit where electrons are "counted")
0069 % 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]
0070 %
0071 % collision mode :
0072 % (0) : Relativistic Maxwellian background
0073 % (1) : High-velocity limit
0074 % (2) : Linearized Belaiev-Budker (momentum-conserving)
0075 %
0076 opt.coll_mode = 0;%2;%
0077 opt.bounce_mode = 0;%1;% bounce-averaged calculation
0078 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0079 %
0080 % normalization mode :
0081 % (0) : "free" mode
0082 % (1) : normalization of f0 maintained to its initial (maxwellian) value. DO NOT use in the runaway problem
0083 % (2) : normalization of f0 maintained to its value at the previous time. Use only in the runaway problem with (+1i)
0084 % (3) : injection at center of particles lost at boundary (use only for stream function)
0085 % (... + 1i) : compensate for numerical error after each time step, acounting for boundary losses.
0086 opt.norm_mode_f = 3;%1i;%1;%0;%3+
0087 %
0088 % time grid (normalized to collision time) -> you can specify :
0089 % - linear arbitrary grid with tn array                 : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0
0090 % - linear arbitrary grid with dtn array                : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10);
0091 % - linear grid with tnmax and grid step                : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000;
0092 % - linear grid with tnmax and number of grid steps     : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i;
0093 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10;
0094 %
0095 % opt.tnmin = 0;% > 0 for log scale
0096 % opt.tn = 10000;% 10000 is time for asymptotic solution
0097 % opt.dtn = 1000;% 10 time steps required for accurate runaway solution - see rundke_dtn
0098 % graph.itdisp = [1,5,10];% display times - only used if opt.timevol > 0
0099 %
0100 % opt.tnmin = 1;% > 0 for log scale
0101 % opt.tn = 1e5;% 10000 is time for asymptotic solution
0102 % opt.dtn = 51;% 10 time steps required for accurate runaway solution - see rundke_dtn
0103 % graph.itdisp = 11:10:51;% display times - only used if opt.timevol > 0
0104 %
0105 %
0106 %
0107 opt.tnmin = 1e0;% > 0 for log scale
0108 %
0109 % opt.tn = [1e8,1e9];%[1e7,1e8];%1e10;% 10000 is time for asymptotic solution
0110 % opt.dtn = [161,90];%[141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn
0111 % graph.itdisp = [61:40:141,251];%[61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0
0112 %
0113 % opt.tn = [1e7,1e8];%1e10;% 10000 is time for asymptotic solution
0114 % opt.dtn = [141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn
0115 % graph.itdisp = [61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0
0116 %
0117 opt.tn = 1e4;%[1e7,1e8];%1e10;% 10000 is time for asymptotic solution
0118 opt.dtn = 81;%[141,90];%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn
0119 graph.itdisp = [21:20:81];%[61:40:141,231];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0
0120 %
0121 %
0122 %
0123 % opt.tnmin = 1;%0;% > 0 for log scale
0124 % opt.tn = 1e8;%1e9;%1e10;% 10000 is time for asymptotic solution
0125 % opt.dtn = 17;%10*1i;%100*1i;%161;%21;%201;% 10 time steps required for accurate runaway solution - see rundke_dtn
0126 % graph.itdisp = [1,10];%[1,10,100];%41:40:161;%5:4:17;%5:4:21;%81:40:201;%1:4:17;% display times - only used if opt.timevol > 0
0127 %
0128 % momentum space grid -> you can specify nmhu_S and either
0129 % - the momentum grid pn_S directly
0130 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows
0131 %
0132 grid.nmhu_S = 141;% number of pitch angle grid points
0133 grid.sfac1 = 80;%40;% % reduction factor for fine grid near abs(mhu) = 1
0134 grid.snmhu1 = 5; % number of fine grid points near abs(mhu) = 1
0135 %
0136 % momentum grid contruction option :
0137 %   (0) do not adjust number of fine and coarse grid points
0138 %   (1) keep fine grid step constant (typically, 0.2)
0139 %       if isreal(np_tail) : keep coarse grid steps constant (typically, 1). Note : do not enforce np_S and np_tail variations
0140 %       if isimag(np_tail), enforce continuity of grid step size at pnmax_S. Note : do not enforce np_S and np_tail variations
0141 %
0142 grid.opt = 1;
0143 %
0144 grid.np_S = 141;% number of bulk momentum grid points
0145 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]
0146 grid.np_tail = 4i;%960;%558;%331;%217;%103;%100;% number of tail momentum grid points
0147 grid.pnmax_S_tail = 1967;%988;%597;%401;205;% maximum of tail momentum (normalized to pT) [107,205,303,401,499,597,793,988,1184,1576,1967] corresponds to [5,10,15,20,25,30,40,50,60,80,100] MeV electrons for betath = 0.1
0148 %
0149 % avalanche modelling :
0150 %
0151 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator
0152 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated
0153 % - For pnmin0_KO = NaN, the limit is set to pnmax_S
0154 %
0155 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator
0156 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy'
0157 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term
0158 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO
0159 %
0160 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator
0161 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used
0162 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference)
0163 %
0164 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator
0165 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons.
0166 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO
0167 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO
0168 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number"
0169 %
0170 opt.avalanche_mode = false;%true;%
0171 opt.pnmin0_KO = NaN;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1
0172 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;%
0173 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator
0174 opt.pnmax2_KO = NaN;%100;%
0175 %
0176 path_simul = '';%Simulation path
0177 %
0178 if opt.bounce_mode == 0,
0179     id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect
0180 else
0181     id_equil = 'TScirc_e1';% toroidal equilibirum with a/R = 1 : then rho_S = r/R
0182 end
0183 id_wave = '';
0184 %
0185 rho_S = 0.5;
0186 %
0187 opt.heavyions = true;
0188 opt.adjustB = true;
0189 %
0190 opt.keyboard = true;
0191 %
0192 graph.tnc = true;
0193 graph.xec = true;
0194 graph.xlog = 1;
0195 graph.xlim = [1e-3,1e2];
0196 graph.xtick = 10.^(-3:2);
0197 %
0198 if alpha == 2,
0199     %
0200     % E/Ec = 2
0201     %
0202     graph.ylim_norm = [1e-25,1e1];
0203     graph.ytick_norm = 10.^(-25:5:5);
0204     graph.ylim_gamma = [1e-26,1e-12];
0205     graph.ytick_gamma = 10.^(-26:4:-12);
0206     graph.ylim_f = [1e-20,1e0];%[1e-45,1e5];
0207     graph.ytick_f = 10.^(-20:4:0);%10.^(-45:5:5);
0208 elseif alpha == 2.5,
0209     %
0210     % E/Ec = 2.5
0211     %
0212     graph.ylim_norm = [1e-15,1e1];
0213     graph.ytick_norm = 10.^(-15:3:0);
0214     graph.ylim_gamma = [1e-25,1e-10];
0215     graph.ytick_gamma = 10.^(-25:3:-10);
0216     graph.ylim_f = [1e-45,1e5];
0217     graph.ytick_f = 10.^(-45:5:5);
0218 elseif alpha == 3,
0219     %
0220     % E/Ec = 3
0221     %
0222     graph.ylim_norm = [1e-11,1e1];%[1e-25,1e1];
0223     graph.ytick_norm = 10.^(-10:2:0);%10.^(-25:5:5);
0224     graph.ylim_gamma = [1e-12,1e-2];
0225     graph.ytick_gamma = 10.^(-12:2:-2);
0226     graph.ylim_f = [1e-20,1e0];%[1e-35,1e5];
0227     graph.ytick_f = 10.^(-20:4:0);%10.^(-35:5:5);
0228     %
0229     graph.dp_cyl = 0;%0.2;
0230     graph.cont = (0:3:100)*1e-11;
0231 %     graph.cont = 1e-7*(-31:2:15)/8;%-3.875e-7:0.25e-7:1.875e-7;%-4e-7:0.25e-7:2e-7;
0232     graph.pnmax_disp = 300;
0233     graph.xlim_disp = [-5,20];
0234     graph.ylim_disp = [0,10];
0235     graph.xtick_disp = -5:5:20;
0236     graph.ytick_disp = 0:5:10;
0237     %
0238 elseif alpha == 5,
0239     %
0240     % E/Ec = 5
0241     %
0242     graph.ylim_norm = [1e-9,1e1];
0243     graph.ytick_norm = 10.^(-10:2:2);
0244     graph.ylim_gamma = [1e-12,1e-4];
0245     graph.ytick_gamma = 10.^(-12:2:-4);
0246 else
0247     graph.ylim_norm = NaN;
0248     graph.ylim_gamma = NaN;
0249 end
0250 %
0251 %************************************************************************************************************************************
0252 %
0253 dp_bulk = grid.pnmax_S/(grid.np_S - 1);
0254 dp_fine = (grid.pnmax_S_tail - grid.pnmax_S)/grid.np_tail;
0255 pfac = imag(grid.np_tail);
0256 if grid.opt == 1,
0257     grid.np_S = 1 + round(grid.pnmax_S/dp_bulk);
0258     if pfac == 0,
0259         grid.np_tail = round((grid.pnmax_S_tail - grid.pnmax_S)/dp_fine);
0260     else
0261         grid.np_tail = round((grid.np_S - 1)*pfac*((grid.pnmax_S_tail/grid.pnmax_S)^(1/pfac) - 1));
0262         grid.opt = pfac;
0263     end
0264 end
0265 %
0266 opt.pnmin_S = sqrt(((1 + Emin*1e3/mc2)^2 - 1)/betath2);
0267 wpr = -taurinv*1i;% real wpr : wpr, positive imaginary part : bhat, negative imaginary part : taurinv.
0268 %
0269 % RUN LUKE
0270 %
0271 [RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,grid,id_equil,id_wave,rho_S);
0272 %
0273 %************************************************************************************************************************************
0274 %
0275 % PROCESS DATA
0276 %
0277 fproc_runaway(RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS,alpha,betath2,wpr,Zi,opt,grid,id_equil,rho_S,graph,p_opt,id_simul);
0278 %
0279 %************************************************************************************************************************************
0280 %
0281 if p_opt > 0,
0282     filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat'];
0283     save(filename);
0284     info_dke_yp(2,['Data saved in ',filename]);
0285 end

Community support and wiki are available on Redmine. Last update: 18-Apr-2019.