rundke_avalanches_tscirc

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