eecoll_dke_yp

LUKE - Electron-electron collision operator

Electron-electron collision operator (1: high velocity limit, 2:
linearized Belaiev-Budker, 3: Lorentz model, 4: non-relativistic
Maxwellian background, 5: Ultrarelativistic Moller model (for runaways
only)

by Yves Peysson (CEA-IRFM,yves.peysson@cea.fr) and Joan Decker (CEA-IRFM,joan.decker@cea.fr) and

0001 function [Fee,Aee,Btee] = eecoll_dke_yp(dkepath,v,z,gamma,sigma,pn,pn2,dpn,J1,J2,spfac,betath_ref,xne_norm,xTe_norm,ir,coll_mode,clustermode)
0002 %LUKE - Electron-electron collision operator
0003 %
0004 %Electron-electron collision operator (1: high velocity limit, 2:
0005 %linearized Belaiev-Budker, 3: Lorentz model, 4: non-relativistic
0006 %Maxwellian background, 5: Ultrarelativistic Moller model (for runaways
0007 %only)
0008 %
0009 %by Yves Peysson (CEA-IRFM,yves.peysson@cea.fr) and Joan Decker (CEA-IRFM,joan.decker@cea.fr) and
0010 %
0011     if nargin < 16,
0012         error('Wrong number of input arguments in eecoll_dke_yp')
0013     end
0014     %
0015     v2 = v.*v;
0016     npn = length(pn);
0017     %
0018     if coll_mode == 1,%High-velocity limit
0019         Fee = xne_norm(ir)'./v2;
0020         Aee = xne_norm(ir)'.*xTe_norm(ir)'./v./v2;
0021         Btee = xne_norm(ir)'./v/2 - xne_norm(ir)'.*xTe_norm(ir)'./v./v2/2;
0022     elseif coll_mode == 0 | coll_mode == 2,%Linearized Belaiev-Budker
0023     %
0024     %Super momentum grid definitions for collision integrals
0025     %
0026     %  spnm = linspace(pn(1)/spfac,pnm(2),spfac);%In order to start from a point not to close to zero
0027     %  for ip=3:npn,
0028     %    spnm = [spnm(1:length(spnm)-1),linspace(pnm(ip-1),pnm(ip),spfac)];
0029     %  end
0030     %
0031 
0032         spn = linspace(pn(1)/spfac,pn(1),spfac);
0033     %  spnp = linspace(pn(1)/spfac,pnp(1),spfac);
0034         for ip=2:npn,
0035             spn = [spn(1:length(spn)-1),linspace(pn(ip-1),pn(ip),spfac)];   
0036     %    spnp = [spnp(1:length(spnp)-1),linspace(pnp(ip-1),pnp(ip),spfac)];
0037         end
0038     %
0039         spn = spn(:);%for numerical integration by trapz_dke_yp subroutine
0040         sdpn = diff([0;spn]);%point i+1/2 (f grid)
0041         spn2 = spn.*spn;
0042         sz = betath_ref*spn;
0043         sz2 = sz.*sz;    
0044         sgamma = sqrt(1 + sz2);
0045         sgamma2 = sgamma.*sgamma;
0046         ssigma = asinh(sz);
0047         sv = spn./sgamma;
0048         sJ1 = -3*sgamma + ssigma.*(3.0./sz + 2*sz);
0049         sJ2 = sgamma - ssigma./sz - 2*sgamma.*sz2/3;
0050      %
0051         sfM = maxwellian_dke_yp(sgamma,spn2,sdpn,xne_norm,xTe_norm,ir);
0052     %
0053         F1_1 = zeros(1,npn);
0054         F1_2 = zeros(1,npn);
0055         F2_1 = zeros(1,npn);
0056     %
0057         B1_1 = zeros(1,npn);
0058         B1_2 = zeros(1,npn);
0059         B1_3 = zeros(1,npn);
0060         B1_4 = zeros(1,npn);
0061         B1_5 = zeros(1,npn);
0062         B2_1 = zeros(1,npn);
0063         B2_2 = zeros(1,npn);
0064         B2_3 = zeros(1,npn); 
0065         %
0066         if isfield(clustermode,'eecoll_dke_yp') && isfield(clustermode.eecoll_dke_yp,'scheduler') && clustermode.eecoll_dke_yp.scheduler.mode,%test distributed computing environment (not usefull for this function)
0067             for ip = 1:npn,arg_range{ip} = ip;end
0068             [dkecluster] = clustermode_luke(clustermode,'eecoll_dke_jd',dkepath);     
0069             [flag,F1_1,F1_2,F2_1,B1_1,B1_2,B1_3,B1_4,B1_5,B2_1,B2_2,B2_3] = mdce_run(@loop_eecoll_dke_yp,{ip,pn,sfM,spn,spn2,sv,sgamma,sgamma2,ssigma,sz,sJ1,sJ2},1,arg_range,dkecluster);
0070             F1_1 = cell2mat(F1_1)';F1_2 = cell2mat(F1_2)';
0071             F2_1 = cell2mat(F2_1)';
0072             B1_1 = cell2mat(B1_1)';B1_2 = cell2mat(B1_2)';B1_3 = cell2mat(B1_3)';B1_4 = cell2mat(B1_4)';B1_5 = cell2mat(B1_5)';
0073             B2_1 = cell2mat(B2_1)';B2_2 = cell2mat(B2_2)';B2_3 = cell2mat(B2_3)';
0074         else
0075             for ip = 1:npn,%Loop for very accurate collision integral calculation close to p = 0 especially (crucial for numerical conservative scheme)
0076                 [F1_1(1,ip),F1_2(1,ip),F2_1(1,ip),B1_1(1,ip),B1_2(1,ip),B1_3(1,ip),B1_4(1,ip),B1_5(1,ip),B2_1(1,ip),B2_2(1,ip),B2_3(1,ip)] = loop_eecoll_dke_yp(ip,pn,sfM,spn,spn2,sv,sgamma,sgamma2,ssigma,sz,sJ1,sJ2);
0077             end
0078         end
0079     %
0080         F1 = F1_1./v2 + F1_2./pn2;
0081         F2 = (1 - sigma./z./gamma).*F2_1./v;
0082         Fee = 4*pi*(F1 + F2);
0083         Aee = 4*pi*xTe_norm(ir)*(F1 + F2)./v;
0084         B1 = B1_1./v/2 - B1_2./v./pn2/6 + B1_3./v./gamma./gamma./z./z/8 - B1_4./v./z./z/4 - B1_5./v./gamma./gamma/4;
0085         B2 = B2_1/2 + B2_2.*(-gamma.*gamma/6 + J1./z./z./gamma/8 - gamma.*J2./z./z/4) - B2_3.*(gamma - sigma./z)./pn2./gamma/4;   
0086         Btee = 4*pi*(B1 + B2);
0087     elseif coll_mode == 3,%Lorentz model
0088         Fee = zeros(1,npn);
0089         Aee = zeros(1,npn);
0090         Btee = zeros(1,npn);
0091     elseif coll_mode == 4,%Maxwellian background
0092         u = v/sqrt(2)/xTe_norm(ir)';
0093         Fee = xne_norm(ir)'.*(erf(u)-2*u.*exp(-u.^2)/sqrt(pi))./v2;
0094         Aee = xne_norm(ir)'.*(erf(u)-2*u.*exp(-u.^2)/sqrt(pi))./v./u./u/2;
0095         Btee = xne_norm(ir)'.*((2*u.^2-1).*erf(u)+2*u.*exp(-u.^2)/sqrt(pi))./v./u./u/4;
0096     elseif coll_mode == 5,%Ultrarelativistic Moller model
0097         Aee = xne_norm(ir)'.*xTe_norm(ir)'./v./v2;
0098         Fee = xne_norm(ir)'./v2;
0099         Btee = xne_norm(ir)'./v/2;
0100     end
0101 %
0102 end
0103 
0104 
0105