frundke_runaway

 
 ***********************This part must be specified by the user, run make files if necessary) *****************************

0001 function [RR,mksa,Ec,f,tn,tRR,tnorm,fgrid,tacc,pmhucrit,ZXXS] = frundke_runaway(id_simul,alpha,betath2,wpr,Zi,opt,varargin)
0002 %
0003 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0004 %
0005 if nargin < 6,
0006    opt = struct;
0007 end
0008 %
0009 if ~isfield(opt,'timevol'),
0010     opt.timevol = 0;
0011 end
0012 %
0013 if ~isfield(opt,'keyboard'),
0014     opt.keyboard = false;
0015 end
0016 %
0017 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants
0018 %
0019 if imag(wpr) ~= 0,%bhat or taurinv provided - ne has to be calculated iteratively because of the coulomb logarithm
0020     %
0021     tol = 1e-6;
0022     nitmax = 20;
0023     if imag(wpr) > 0,% input wpr is bhat
0024         bhat2 = imag(wpr)^2;
0025     elseif imag(wpr) < 0,% wpr is taurinv
0026         bhat2 = -imag(wpr)*betath2^(3/2);
0027     end
0028     %
0029     [~,equil,dkeparam] = fprep_runaway(0,Zi,opt,varargin{:});% magnetic field and ne20 first guess
0030     pB = sqrt(equil.ptBx(:,1).^2 + equil.ptBy(:,1).^2 + equil.ptBPHI(:,1).^2).';
0031     pne20 = equil.pne/1e20;
0032     prho = equil.ptx(:,1)/equil.ptx(end,1);
0033     %
0034     B = interp1(prho,pB,dkeparam.rho_S);
0035     ne20 = interp1(prho,pne20,dkeparam.rho_S);
0036     %
0037     for it = 1:nitmax,
0038         lnc_e_ref = 31.3 - 0.5*log(ne20*1e20) + log(mc2*1000) + log(betath2);%Reference Coulomb logarithm (Sauter et al. Phys. Plasmas, 6 (1999) 2834)
0039 %         lnc_e_ref = 17.6825;
0040         wpr = bhat2*(3/2)*lnc_e_ref/betath2^(3/2);
0041         %
0042         ne20_old = ne20;
0043         ne20 = (B*qe/me)^2*me*e0/(qe^2*wpr)/1e20;
0044         %
0045         convpar = abs(ne20 - ne20_old)/ne20_old;
0046         %
0047         if (isfield(opt,'adjustB') && opt.adjustB) || convpar <= tol,
0048             break
0049         end
0050     end
0051     %
0052     if convpar > tol,
0053         disp('WARNING : no convergence found for ne calculation from bhat')
0054     end
0055 end
0056 %
0057 [dkepath,equil,dkeparam,dkedisplay,waves] = fprep_runaway(wpr,Zi,opt,varargin{:});
0058 %
0059 dkeparam.timevol = opt.timevol;% time evolution (1) or steady-state solution (0)
0060 dkeparam.keyboard = opt.keyboard;% give hand if negative density (1) or return (0)
0061 %
0062 equil.pTe = betath2*mc2*ones(size(equil.pTe));
0063 %
0064 epsi = alpha*betath2;
0065 ohm = ohm_flat(equil,epsi);
0066 %
0067 [Znorm,~,~,dke_out,~,~,momentumDKE,~,Zmomcoef,~,~,mksa] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,'','',[],[]);
0068 %
0069 Ec = Zmomcoef.Ec_ref/1e3;% energy in MeV
0070 %
0071 fgrid = selectfields_jd(momentumDKE,{'pn','pnp','pn_S','dpn','mhu','mhu_S','dmhu'});
0072 ZXXS = dke_out.ZXXS;
0073 %
0074 if opt.timevol == 0,
0075     %
0076     if dkeparam.norm_mode_f == 0,
0077         RR = dke_out.xRRp_fsav/Znorm.x_0_fsav;
0078     else
0079         RR = dke_out.xRRp_fsav;
0080     end
0081     %
0082     f = [dke_out.XXfM,dke_out.XXf0(end)];
0083     %
0084     tn = dke_out.FPtime(end);
0085     tnorm = '';
0086     tRR = '';
0087     %
0088     tacc = dke_out.XXRintmask(end);
0089     %
0090 else
0091     %
0092     tn = [0,dke_out.FPtime];
0093     dtn = diff(tn);
0094     ntn = length(dtn);
0095     %
0096     if isempty(dke_out.xRRp_fsav{end}),% the calculation stopped before the end (negative norm at some point)
0097         RR = NaN;
0098     elseif dkeparam.norm_mode_f == 0,
0099         RR = dke_out.xRRp_fsav{end}/Znorm(end).x_0_fsav;
0100     else
0101         RR = dke_out.xRRp_fsav{end};
0102     end
0103     %
0104     f = [dke_out.XXfM,dke_out.XXf0];
0105     %
0106     tnorm.re = zeros(1,ntn+1);
0107     tnorm.ri = NaN(1,ntn+1);
0108     tnorm.i = ones(1,ntn+1);
0109     %
0110     tacc = dke_out.XXRintmask;
0111     %
0112     for itn = 1:ntn,
0113         if length(Znorm) >= itn,% the calculation ran to this point
0114             tnorm.re(itn+1) = dke_out.xnorm_rem{itn} + dke_out.xnorm_rep{itn};
0115             tnorm.ri(itn+1) = Znorm(itn).x_Rfrac_0_fsav;
0116             tnorm.i(itn+1) = Znorm(itn).x_0_fsav;
0117         else% the calculation stopped before the end (negative norm at some point)
0118             tnorm.re(itn+1) = NaN;
0119             tnorm.ri(itn+1) = NaN;
0120             tnorm.i(itn+1) = NaN;
0121             %
0122             f{itn+1} = NaN(size(f{1}));
0123             %
0124             tacc{itn} = NaN(size(f{1}));
0125         end            
0126     end
0127     %
0128     tnorm.b = tnorm.i - tnorm.ri;
0129     tnorm.tot = tnorm.i + tnorm.re;% to check particle conservation
0130     %
0131     tRR.re = diff(tnorm.re)./dtn./tnorm.b(2:end);
0132     tRR.ri = diff(tnorm.ri)./dtn./tnorm.b(2:end);
0133     %
0134 end
0135 %
0136 XXFp_ij = dke_out.ZXXF_c.p_ij + dke_out.ZXXF_e.p_ij + dke_out.ZXXF_s.p_ij;
0137 np = length(fgrid.pn);
0138 pmhucrit = NaN(1,np);
0139 for ipn = 1:np,
0140     imhucrit = find(XXFp_ij(ipn,1:end-1).*XXFp_ij(ipn,2:end) < 0);
0141     if isempty(imhucrit),
0142         if XXFp_ij(ipn,end) > 0,
0143             pmhucrit(ipn) = -1;
0144         elseif XXFp_ij(ipn,end) < 0,
0145             pmhucrit(ipn) = 1;
0146         end
0147     elseif length(imhucrit) == 1,% linear fit
0148         pmhucrit(ipn) = fgrid.mhu(imhucrit) + (fgrid.mhu(imhucrit + 1) - fgrid.mhu(imhucrit))*XXFp_ij(ipn,imhucrit)/(XXFp_ij(ipn,imhucrit) - XXFp_ij(ipn,imhucrit + 1));
0149     end
0150 end
0151 
0152