runaway_sinksource

 Calculates the Runaway Secondary (Avalanche) Source and Loss Terms, bounce-averaged or not
 the results are NOT normalized.         15/06/2007  modified  11/07/07

0001 function [XXRS,xRL]=runaway_sinksource(dkeparam,mksa,gridDKE,equilDKE,Zbouncecoef,Zmomcoef,xne,xnr_loss,ohm,Xxpnlim)
0002 
0003 % Calculates the Runaway Secondary (Avalanche) Source and Loss Terms, bounce-averaged or not
0004 % the results are NOT normalized.         15/06/2007  modified  11/07/07
0005 
0006 % INPUT
0007 %
0008 %     - dkeparam : simulation parameters
0009 %     - mksa : normalization constants
0010 %     - gridDKE : grid parameters
0011 %     - equilDKE : magnetic equilibrium interpolated on the resolution grid
0012 %     - Zbouncecoef : used to get lambda and qtilde, (bounce-averaging parameters)
0013 %     - Zmomcoef : used to get gamma and XXgamma
0014 %     - ohm : electric field (used to know the propagation direction of the runaway (opposed to the
0015 %                            electric field)
0016 %     - xnr_loss : density of already existing runaway electrons (1 nr)
0017 %     - xne : density of thermal electron likely to be collided (inside Xxpnlim)
0018 %     - Xxpnlim : momentum separating the source and loss domain of the
0019 %                collision process (purely conventionnal) (nmhu nr)
0020 %
0021 
0022 % OUTPUT
0023 %
0024 %     - XXSR : Runaway Secondary Source  (npn, nmhu, nr)
0025 %     - xRL  : Runaway Secondary Loss of the termal bulk  (1,nr)
0026 %
0027 % the theta grid has to be uniform (if not, a line must be changed)
0028 %
0029 % The loss term is calculated as the integral of the source term on the
0030 % outward domain delimited by Xxpnlin
0031 %
0032 % The source term is given everywhere. Caution : the bounce-averaged sounce
0033 % term is not rigourously correct : the symmetry of the trapped region has
0034 % still to be enforced (cf the note on the report). This has been done
0035 % in sinksourcecoeffbuilder.m .
0036 %
0037 if nargin < 8,
0038     error('Not enough input arguments');
0039 end
0040 %
0041 %    Initialisation
0042 %
0043 bounce_mode=dkeparam.bounce_mode;
0044 %
0045 betath_ref=mksa.betath_ref;
0046 lnc_e_ref=mksa.lnc_e_ref;
0047 nhu_ref=mksa.nhu_ref;
0048 ne_ref=mksa.ne_ref;
0049 %
0050 XXgamma=Zmomcoef.XXgamma;
0051 XXpn=gridDKE.XXpn;
0052 XXmhu=gridDKE.XXmhu;
0053 gamma=Zmomcoef.gamma;
0054 mhu=gridDKE.mhu;
0055 npn=gridDKE.npn;
0056 nmhu=gridDKE.nmhu;
0057 nr=gridDKE.nr;
0058 %
0059 XXnebar=permute(repmat(xne/mksa.ne_ref,[npn 1 nmhu]),[1 3 2]);
0060 XXnrbar=permute(repmat(xnr_loss/mksa.ne_ref,[npn 1 nmhu]),[1 3 2]);
0061 %
0062 dir=ohm.epsi(1)/abs(ohm.epsi(1)); %runaway propagation direction  CAUTION
0063 % it depends on the convention on the sign of the E field. Here, runaway electrons
0064 % moves towards mhu=+1 for E > 0.
0065 %
0066 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0067 % bounce-averaged source term
0068 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0069 if bounce_mode==1   
0070     %
0071     XXthetaT=zeros(npn,nmhu,nr);
0072     XXRS=zeros(npn,nmhu,nr);
0073     %
0074     XXB0=permute(repmat(equilDKE.xB0,[npn 1 nmhu]),[1 3 2]);
0075     ap=equilDKE.ap;
0076     XXqtilde=permute(repmat(Zbouncecoef.xqtilde,[npn 1 nmhu]),[1 3 2]);
0077     XXlambda=permute(repmat(Zbouncecoef.Xxlambda,[1 1 npn]),[3 1 2]);  %Xxlambda was (nmhu,nr)
0078     %
0079     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0080     % Theta star depending term calculation
0081     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0082     %
0083     theta_S=equilDKE.mtheta;
0084     ntheta_S=length(theta_S);
0085     %
0086     dtheta=zeros(1,ntheta_S-1);
0087     dtheta=theta_S(2:end)-theta_S(1:end-1);
0088     %
0089     XXBstar=2*XXB0./( (1-XXmhu.^2).*(1+XXgamma)); %field to be interpolated
0090     YxBstar=reshape(XXBstar,[npn*nmhu nr]);
0091     %
0092     XB=sqrt(equilDKE.XBx.^2+equilDKE.XBy.^2+equilDKE.XBphi.^2);
0093     %
0094     xBmax=XB(:,(ntheta_S+1)/2); %Attention (nr,1)
0095     YxBmax=permute(repmat(xBmax,[1 npn*nmhu]),[2 1]); %(np*nmhu nr)
0096     YxB0=reshape(XXB0,[npn*nmhu nr]);
0097     %
0098     Ymhu=reshape(XXmhu,[1,npn*nmhu]); %for direction selection
0099     %
0100     for ir=1:nr
0101         %
0102         YthetaT=zeros(1,npn*nmhu); %to be used as transfer variable in the interpolation loop
0103         %
0104         Ymask=find(YxBmax(:,ir)>=YxBstar(:,ir) & YxBstar(:,ir)>=YxB0(:,ir) & Ymhu(:)*dir>=0 ); % Existence domain of thetastar for a given value of r
0105         %
0106         YB=YxBstar(Ymask,ir);
0107         Bx=equilDKE.XBx(ir,:);
0108         By=equilDKE.XBy(ir,:);
0109         Bphi=equilDKE.XBphi(ir,:);  
0110         x=equilDKE.Xx(ir,:);
0111         y=equilDKE.Xy(ir,:);
0112         %
0113         %
0114         % Matrix to be interpolated at thetastar
0115         r_S=sqrt(x.^2+y.^2);   % (1;ntheta_S)
0116         Bp_S=sqrt(Bx.^2+By.^2);  % (1;ntheta_S)
0117         Psipr_S=abs(Bx.*y - By.*x)./(Bp_S.*r_S);  % (1;ntheta_S)
0118         %
0119         B0=equilDKE.xB0(ir); % scalar
0120         Psi_S=sqrt(Bx.^2+By.^2+Bphi.^2)./B0;% (1;ntheta_S)
0121         dPsidtheta_S=zeros(1,ntheta_S);
0122         dPsidtheta=zeros(1,ntheta_S);
0123         dPsidtheta(2:end)=(Psi_S(2:end)-Psi_S(1:end-1))./(dtheta); %on the half grid
0124         dPsidtheta(1)=dPsidtheta(end); % periodicity
0125         dPsidtheta_S(1:end-1)=(dPsidtheta(1:end-1)+dPsidtheta(2:end))/2; 
0126         % here, the full grid is in the middle of the half grid
0127         % as the theta grid is uniform
0128         dPsidtheta_S(end)=dPsidtheta_S(1); %periodicity
0129         %
0130         YT=zeros(ntheta_S,1,4);  %the second 1 is nr, here =1 beacause of the loop - to be adjusted
0131         YT(:,1,1)=r_S(:);
0132         YT(:,1,2)=Bp_S(:);
0133         YT(:,1,3)=Psipr_S(:);
0134         YT(:,1,4)=dPsidtheta_S(:); %(4, ntheta)
0135         %
0136         % Each term, and especially dPsidtheta, has to be interpolated
0137         % separately in order to preserve the source-term integrability.
0138         % Do not interpolate directly 1/dPsidtheta
0139         %
0140         [Ytheta_star1,Ytheta_star2,YT1,YT2] =...
0141         interp_theta_jd(theta_S,Bx',By',Bphi',YB',YT);    
0142         %
0143         r1=YT1(1,:,1);
0144         Bp1=YT1(1,:,2);
0145         Psipr1=YT1(1,:,3);
0146         dPsidtheta1=YT1(1,:,4);
0147         %
0148         r2=YT2(1,:,1);
0149         Bp2=YT2(1,:,2);
0150         Psipr2=YT2(1,:,3);
0151         dPsidtheta2=YT2(1,:,4);
0152         %
0153         YmaskthetaT1=r1./(Bp1.*Psipr1.*abs(dPsidtheta1)); 
0154         YmaskthetaT2=r2./(Bp2.*Psipr2.*abs(dPsidtheta2)); 
0155         %
0156         YmaskthetaT = YmaskthetaT1 + YmaskthetaT2; 
0157         YthetaT(Ymask)=YmaskthetaT;
0158         XXthetaT(:,:,ir)= reshape(YthetaT,[npn nmhu]);
0159         %
0160     end
0161     %
0162     % Preliminary Constant Calcutation
0163     C = (nhu_ref*ne_ref)*(1/(2*pi))*(1/(pi*lnc_e_ref*ap));
0164     %
0165     % Final bounce-averaged result
0166     XXRS=C*(XXB0./(XXlambda.*XXqtilde)).*XXnebar.*XXnrbar.*(1./(XXgamma.*XXpn.^3.*(XXgamma-1))).*( abs(XXmhu)./((1-XXmhu.^2).^2)) .*XXthetaT;
0167     %
0168 else
0169     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0170     % source term without bounce-average
0171     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0172     %
0173     XXmhustarval= dir * sqrt((XXgamma-1)./(XXgamma+1));
0174     XXmhustar=interp1(mhu,mhu,XXmhustarval,'nearest');    
0175     XXdeltastar=(XXmhu==XXmhustar);        
0176     XXdmhu=repmat(gridDKE.dmhu,[npn 1 nr]);
0177     %
0178     XXRS=(nhu_ref*ne_ref)*(betath_ref^2)/(4*pi*lnc_e_ref).*XXnebar.*XXnrbar.*(1./(XXgamma.*XXpn.*(XXgamma-1).^2)).*XXdeltastar./XXdmhu;
0179     %
0180 end
0181 
0182 %%%%%%%%%% Loss term, calculated using partially analytical integration %%%%%%%%%
0183 %
0184 xRL=zeros(1,nr);
0185 mhustarval= dir*sqrt( (gamma-1)./(gamma+1) );
0186 mhustar=interp1(mhu(1:end),[1:nmhu]',mhustarval,'nearest');
0187 %
0188 %
0189 for ir=1:gridDKE.nr
0190     H=(gridDKE.pn(:)>=Xxpnlim(mhustar(:),ir));
0191     if sum(H)>0
0192         xRL(ir)=(nhu_ref*ne_ref)*((betath_ref^(5)*(xnr_loss(ir)/ne_ref)*(xne(ir)/ne_ref))/(2*lnc_e_ref))*integral_dke_jd(gridDKE.dpn(H),(gridDKE.pn(H))./(gamma(H).*(gamma(H)-1).^2),2); 
0193     else
0194         xRL(ir)=0;
0195     end
0196 end
0197 
0198 
0199 %%%%% Complementary loss term %%%%%%%%%%%%%%%%%%%%%
0200 % Loss term
0201 xRLC=zeros(1,nr);
0202 xRLC(:)=(nhu_ref*(mksa.betath_ref)^3*xne(:).*xnr_loss(:)./(2*mksa.lnc_e_ref*ne_ref)).*(1./(sqrt(1+(mksa.betath_ref*dkeparam.pnmax_S)^2)-1) ) ;
0203 
0204 xRL=xRL+xRLC;
0205 
0206 %if sum(xRL<=0)>0 keyboard; end