rundke_avalanches

 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 for Zii = 1,
0012 %
0013 permission = test_permissions_yp;
0014 %
0015 if ~permission 
0016     disp('Please move the script to a local folder where you have write permission before to run it')
0017     return;
0018 end
0019 %
0020 % grid parameters
0021 %
0022 pnmax_S = 28;%20;% for E/Ec = 4, Pc = 5; for E/Ec = 6, Pc = 4;
0023 np_S = 201;%201;%
0024 nmhu_S = 101;
0025 f_pmax = floor(100/pnmax_S); %1;% %for avalanche calculation. if avalanche_mode and not specified = 1
0026 %
0027 mfactor = 100;%Factor for super theta grid (equil -> equilDKE)
0028 %
0029 %
0030 % normalized electric field (with respect to the Critical field Ec)
0031 %
0032 % temperature betath2 = Te/mc2
0033 % (betath2 = 10^(-6) is validated for NR limit)
0034 %
0035 betath2 = 1e-2;%1e-5;%
0036 %
0037 alpha = 4; %E/Ec
0038 epsi = alpha*betath2;
0039 %
0040 pnmin_S = NaN;%sqrt(2)/sqrt((alpha-1)*betath2);%5;%1/sqrt((alpha-1)*betath2);%NaN;%
0041 %
0042 % ion charge (single species, cold ions)
0043 %
0044 Zi = Zii;%
0045 %
0046 % collision mode :
0047 % (0) : Relativistic Maxwellian background
0048 % (1) : High-velocity limit
0049 % (2) : Linearized Belaiev-Budker (momentum-conserving)
0050 %
0051 coll_mode = 0;
0052 %
0053 id_simul = 'Runaway';%Simulation ID
0054 %
0055 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0056 %
0057 %id_simul = 'Runaway';%Simulation ID
0058 %path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0059 %
0060 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0061 rho_S = [0.1];%1%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0062 %
0063 id_path = '';%For all paths used by DKE solver
0064 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0065 %
0066 id_equil = 'TScyl';%For plasma equilibrium
0067 %id_equil = 'TScirc';%For plasma equilibrium
0068 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0069 %
0070 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0071 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0072 %
0073 id_display = 'PARTIAL_VISUAL';%'NO_DISPLAY''PARTIAL_VISUAL';%For output code display, files from Project_DKE/DISPLAY_FILES/
0074 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0075 %
0076 id_transpfaste = 'magneticturbulence';
0077 path_transpfaste = '';
0078 %************************************************************************************************************************************
0079 %************************************************************************************************************************************
0080 %************************************************************************************************************************************
0081 %
0082 [dkepath,equil,dkeparam,dkedisplay,dummy,dummy,transpfaste] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'','','','','','','transpfaste',id_transpfaste,path_transpfaste);
0083 %
0084 %**************************************************************************
0085 %**************************************************************************
0086 %
0087 if exist('dmumpsmex','file') == 3,
0088     dkeparam.invproc = -2;
0089 end
0090 %
0091 dkeparam.pnmin_S = pnmin_S;
0092 dkeparam.pnmin2_KO = 1*1i; % low threshold on secondary electrons for the knock-on operator
0093 dkeparam.pnmax2_KO = pnmax_S*f_pmax; % high threshold on secondary electrons for the knock-on operator (>= pnmax_S)
0094 dkeparam.pnmax_S = pnmax_S;
0095 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0096 dkeparam.norm_mode_f = 0;% + 1i;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0097 dkeparam.timevol = 1;
0098 %dkeparam.nit_ohm = 2;
0099 dkeparam.tn = 30000;%time for asymptotic solution with norm_mode_f = 1
0100 dkeparam.dtn = 1000;% time steps required for accurate runaway solution - see rundke_dtn
0101 dkeparam.coll_mode = coll_mode;% mode for the collision operator
0102 dkeparam.bounce_mode = 0; %bounce average
0103 dkeparam.avalanche_mode = 1; %if 1 avalanche mode on, 0 or undefined for off
0104 %
0105 options_fzero = optimset('fzero');
0106 options_fzero.Display = 'iter';
0107 options_fzero.TolX = 1e-10;%convergence tolerance
0108 dkeparam.options_fzero = options_fzero;
0109 %
0110 transpfaste.Dr0 = 0;
0111 %transpfaste.vparmin = 1;
0112 %
0113 equil.pzTi = 1e-10*ones(size(equil.pzTi));% cold ions
0114 %
0115 [qe,me,dummy,~,e0,dummy,dummy,mc2,clum] = pc_dke_yp;%Universal physics constants
0116 %
0117 equil.pTe = betath2*mc2*ones(size(equil.pTe));
0118 %
0119 equil.zZi = [1,1,1,Zi];
0120 equil.pzni = [0;0;0;1/Zi]*equil.pne;% single ion species
0121 %
0122 epsi = alpha*betath2;
0123 %
0124 xnre_init = 0;%zeros(1,length(rho_S)); const. value, should be changed to vectorial form(one for each rho)
0125 %
0126 ohm = ohm_flat(equil,epsi);
0127 %
0128 ohm.xnre_init = xnre_init;
0129 %
0130 dkeparam.pn_S = linspace(0,pnmax_S,np_S);
0131 dkeparam.nmhu_S = nmhu_S;
0132 dkeparam.psin_S = psin_S;
0133 dkeparam.rho_S = rho_S;
0134 dkeparam.mfactor = mfactor;
0135 %
0136 % Global parameters for the vectorial magnetic equilibrium
0137 %
0138 fitparam.mode_equil = 1;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y)
0139 fitparam.nharm = NaN;%Number of harmonics in the magnetic equilibrium interpolation (less than ntheta_equil/2)
0140 fitparam.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters
0141 fitparam.method = 'spline';%'linear';%nearest,spline,pchip
0142 %
0143 equil.equil_fit = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);%Build vectorized magnetic equilibrium structure
0144 %
0145 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,'',transpfaste,'',[],[]);
0146 %
0147 RR_an = RR_rosenbluth_en(alpha, Zi, mksa.ne_ref, mksa.lnc_e_ref, equilDKE.ap*equilDKE.xrho/equilDKE.Rp); %analytic AVALANCHE growth rate Rosenbluth Putvinski NF 1997
0148 %
0149 %RR = dke_out.XxRR_fsav(end,:);
0150 Rfrac = Znorm(end).x_Rfrac_0_fsav; %runaway fraction (nr of runaways between pmin and pmax) primary + avalanche
0151 %
0152 xnre = cell2mat(dke_out.xnorm_rep);
0153 %
0154 RRAint = Znorm(1,end).xRRAint./xnre(end-1); %integral over S[pmin pmax]
0155 RRAinttot = Znorm(1,end).xRRAtot./xnre(end-1); %total integral S[0 pmax]
0156 RR = dke_out.XxRR_fsav{end}; %runaway loss rate (primary)
0157 RR = RR(end,:); %[vol-1 tau_c^-1]
0158 %RRunnorm = RR*mksa.nhu_ref; %vol-1 s^-1
0159 %RRRP = RR*mksa.nhu_ref/nre(1,1);
0160 for it = 1:(length(xnre)-1)
0161     gr(it) = (xnre(it+1) - xnre(it))/xnre(it)/dkeparam.dtn;
0162 end
0163 figure;graph1D_jd(1:(length(xnre)-1),gr,0,0,'iteration','growth rate (s^-1)','Numerical growth rate',NaN,NaN,NaN,'-','o','b',1)
0164 hold on
0165 %plot(1:(length(xnre)-1), RRAinttot,'rx')
0166 %plot(1:(length(xnre)-1), RRAint,'kx'); hold off
0167 %
0168 if dkeparam.avalanche_mode == 1
0169     pnext = pnmax_S*f_pmax;
0170     gamma_max = sqrt((pnmax_S*mksa.betath_ref)^2+1);
0171     gamma_ext = sqrt((pnext*mksa.betath_ref)^2+1);
0172     %Ec = mksa.ne_ref*qe^3*mksa.lnc_e_ref/(4*pi*e0^2*me*clum^2);
0173     %E = alpha * Ec;
0174     %nhu_c = qe^4*mksa.ne_ref*mksa.lnc_e_ref/(4*pi*e0^2*me^2*clum^3);
0175     %nhu_c = mksa.nhu_ref*mksa.betath_ref^3;
0176     dtn_min = me*clum*(sqrt(gamma_ext^2-1)-sqrt(gamma_max^2-1))/(qe*mksa.Edreicer_ref*epsi)*mksa.nhu_ref; %mininmum time step dtn for dreicer acc from pmax to pext
0177    
0178     %dtn_min = me*clum*(sqrt(gamma_ext^2-1)-sqrt(gamma_max^2-1))/(qe*E)*mksa.nhu_ref; %mininmum time step dtn for dreicer acc from rest to gamma
0179     gam = gamma_max:0.01:(gamma_max+10);
0180     dt_acc = me*clum*(sqrt(gam.^2-1)-sqrt(gamma_max^2-1))/(qe*mksa.Edreicer_ref*epsi)*mksa.nhu_ref; %mininmum time step dtn for dreicer acc from rest to gamma
0181     titl=['(ro): gamma @ p_{max} = ',num2str(pnmax_S),'p_{th}'];
0182     figure;graph1D_jd(gam,dt_acc,0,0,'gamma','dt_n Dreicer',titl,NaN,NaN,NaN,'-','none','b',2)
0183     hold on
0184     graph1D_jd(gamma_ext,dtn_min,0,0,'','','',NaN, NaN,NaN,'none','o','r',3);
0185     hold off
0186     disp(['     Mininmum time step dtn for dreicer acc from rest to gamma(p_max) = : ',num2str(gamma_max),' is tau*nhu_ref = ',num2str(dtn_min)])
0187 end
0188 %
0189 if length(rho_S) == 1
0190 %[logXf0_cyl_ref,ppar_cyl_ref,dppar_cyl_ref,pperp_cyl_ref,dpperp_cyl_ref] = s2c_dke_yp(log(cell2mat(dke_out.XXf0)),momentumDKE.pn,momentumDKE.mhu,0.1);%Spherical to cylindrical coordinate transformation
0191 [logXf0_cyl_ref,ppar_cyl_ref,dppar_cyl_ref,pperp_cyl_ref,dpperp_cyl_ref] = s2c_dke_yp(log(dke_out.XXf0{end}),momentumDKE.pn,momentumDKE.mhu,0.1);%Spherical to cylindrical coordinate transformation
0192 Xf0_cyl_ref = exp(logXf0_cyl_ref);%For accurate representation in figures
0193 %
0194 logXfinit_cyl_ref = s2c_dke_yp(log(dke_out.XXfinit),momentumDKE.pn,momentumDKE.mhu,0.1);%Spherical to cylindrical coordinate transformation
0195 Xfinit_cyl_ref = exp(logXfinit_cyl_ref);%For accurate representation in figures
0196 %
0197 Xppar_cyl_ref = ones(length(pperp_cyl_ref),1)*ppar_cyl_ref(:)';%Grid for the distribution functions only
0198 Xpperp_cyl_ref = pperp_cyl_ref(:)*ones(1,length(ppar_cyl_ref));
0199 Xp_cyl_ref = sqrt(Xppar_cyl_ref.*Xppar_cyl_ref + Xpperp_cyl_ref.*Xpperp_cyl_ref);
0200 %
0201 Xf0_cyl_ref = Xf0_cyl_ref.*(Xf0_cyl_ref>0);%Remove negative values
0202 Xf0_cyl_ref = Xf0_cyl_ref.*(Xp_cyl_ref<=max(ppar_cyl_ref));%Remove values above max(ppar)
0203 Xf0_cyl_ref(isnan(Xf0_cyl_ref)) = 0;%Remove NaN values
0204 %
0205 Xfinit_cyl_ref = Xfinit_cyl_ref.*(Xfinit_cyl_ref>0);%Remove negative values
0206 Xfinit_cyl_ref = Xfinit_cyl_ref.*(Xp_cyl_ref<=max(ppar_cyl_ref));%Remove values above max(ppar)
0207 Xfinit_cyl_ref(isnan(Xfinit_cyl_ref)) = 0;%Remove NaN values
0208 %
0209 % Remove values less than computer accuracy
0210 %
0211 Xmask = abs(Xf0_cyl_ref) < eps;
0212 %
0213 pnmax = max(momentumDKE.pn_S);
0214 pnmax_0 = min(min(Xp_cyl_ref(Xmask)));
0215 %
0216 Xmask = Xp_cyl_ref >= pnmax;    
0217 %
0218 Xf0_cyl_ref(Xmask) = 0;
0219 Xfinit_cyl_ref(Xmask) = 0;
0220 %
0221 figure(1),clf
0222 %
0223 x1 = ppar_cyl_ref;
0224 xlab1 = 'p_{||}/p_{Te}';
0225 xlim1 = [-pnmax,pnmax];
0226 %
0227 x2 = pperp_cyl_ref;
0228 xlab2 = 'p_{\perp}/p_{Te}';
0229 xlim2 = [0,pnmax];
0230 %
0231 leg = {'f_M','f'};
0232 %
0233 y1 = Xfinit_cyl_ref;
0234 y2 = Xf0_cyl_ref;
0235 %
0236 y9 = sqrt(pnmax^2 - x1.^2);
0237 %
0238 figure;
0239 graph2D_jd(x1,x2,y1','','','',xlim1,xlim2,0,mksa.betath_ref,[0:0.5:pnmax],0,'-','b',0.5,20+14i,1,max(max(y1)));
0240 hold on
0241 graph2D_jd(x1,x2,y2','','','',xlim1,xlim2,0,mksa.betath_ref,[0:0.5:pnmax],0,'-','r',2,20+14i,1,max(max(y1)));
0242 hold off
0243 %
0244 graph1D_jd(x1,y9,0,0,xlab1,xlab2,'',NaN,xlim1,xlim2,'-','none','k',2,20+14i,gca,0.9,0.7,0.7);
0245 %
0246 axis('equal');axis([xlim1,xlim2]);
0247 legend(leg);
0248 end
0249 %
0250 disp('Simulation parameters : ')
0251 disp(['coll_mode = ',num2str(coll_mode),' ; Te/mc2 = ',num2str(betath2),' ; E/Ec = ',num2str(alpha),' ; E/ED = ',num2str(alpha*betath2),' ; Zi = ',num2str(Zi)])
0252 disp(['NR limit : pc/pT = ',num2str(1/sqrt(alpha*betath2)),' ; pmax/pc = ',num2str(pnmax_S*sqrt(alpha*betath2))])%NR
0253 disp(['Relativ. : pc/pT = ',num2str(1/sqrt((alpha - 1)*betath2)),' ; pmax/pc = ',num2str(pnmax_S*sqrt((alpha - 1)*betath2))])%rel
0254 disp(' ')
0255 disp(['--> LUKE RR (normalized to ne and tauc) : ',num2str(RR)])
0256 disp('----------------------')
0257 %disp(['--> RR (normalized to ne) : ',num2str(RRunnorm)])
0258 %disp('----------------------')
0259 %
0260 diary off
0261 %
0262 %figure;graph1D_jd((0:dkeparam.dtn:dkeparam.tn),xnre,0,0,'t_n','number of runaways (norm to n_e)','',NaN,NaN,NaN,'-','none','b',2)
0263 %figure;graph1D_jd((0:dkeparam.dtn:dkeparam.tn),xnre,0,1,'t_n','number of runaways (norm to n_e)','',NaN,NaN,NaN,'-','none','b',2)
0264 %
0265 %curr0 = cell2mat(dke_out.curr0'); %*mksa.j_ref
0266 %figure;graph1D_jd((1:dkeparam.dtn:dkeparam.tn),curr0,0,0,'t_n','J_0','',NaN,NaN,NaN,'-','none','b',2)
0267 %figure;graph1D_jd((1:dkeparam.dtn:dkeparam.tn),cell2mat(dke_out.xcurr),0,0,'t_n','J_{tot}','',NaN,NaN,NaN,'-','none','b',2)
0268 %
0269 %growthrate(ia) = gr(end)*mksa.nhu_ref
0270 %
0271 figure
0272 graph1D_jd((1:dkeparam.dtn:dkeparam.tn),xnre(:),0,0,'t_n','n_{re}/n_{tot}','',NaN,NaN,NaN,'-','none','b',2)
0273 figure
0274 grDA = diff(xnre(:))'./dkeparam.dtn./(1-xnre(1:(end-1)));
0275 graph1D_jd(xnre(1:(end-1)),grDA,0,0,'n_r','1/(1-n_r)dn_r/dt_n','',NaN,NaN,NaN,'-','none','b',2)
0276 grA = (grDA(end)-grDA(end-20))/(xnre(end)-xnre(end-20));
0277 %
0278 %************************************************************************************************************************************
0279 %
0280 %
0281 %
0282 %eval(['save ','DKE_RESULTS_',id_equil,'_D_beta2_5e-3_EEc_',num2str(alpha),'.mat']);
0283 %info_dke_yp(2,['Data saved in ','DKE_RESULTS_',id_equil,'_DA','_rho_0',num2str(irho),'_tn_',num2str(dkeparam.tn),'_E_',num2str(alpha),'.mat']);
0284 %eval(['save','PROC_',id_equil,'_DA_','beta2_1e-3','_EEc_',num2str(alpha),'xnre','alpha','dke_out','mksa']);
0285 %save(['PROC_TScyl_DA_beta2_2e-5_EEc_',num2str(alpha)],'xnre','alpha','dke_out','mksa')
0286 %eval(save('PROC_TScyl_DA_beta2_1e-2_EEc_10','xnre','alpha','dke_out','mksa'))
0287 %info_dke_yp(2,['Data saved in ','DKE_RESULTS_',id_equil,'_DA_','beta2_1e-3','_EEc_',num2str(alpha),'.mat']);
0288 %eval(['save ','E_50_Z',num2str(Zi),'.mat']);
0289 %info_dke_yp(2,['Data saved in',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'_',filename,'.mat']);
0290 end