rad_dke_yp

PURPOSE ^

SYNOPSIS ^

function [ZXXD_a,ZXXF_a] = rad_dke_yp(dkeparam,display,equilDKE,gridDKE,mksa,Zbouncecoef,Zmomcoef,transpfaste)

DESCRIPTION ^

 Calculation of the anomalous radial fast electron transport matrix coefficients for arbitrary magnetic plasma equilibrium 

 by Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr> and J. Decker (CEA-IRFM) <joan.decker@cea.fr>

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ZXXD_a,ZXXF_a] = rad_dke_yp(dkeparam,display,equilDKE,gridDKE,mksa,Zbouncecoef,Zmomcoef,transpfaste)
0002 %
0003 % Calculation of the anomalous radial fast electron transport matrix coefficients for arbitrary magnetic plasma equilibrium
0004 %
0005 % by Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr> and J. Decker (CEA-IRFM) <joan.decker@cea.fr>
0006 %
0007 if nargin < 7,
0008     error('Wrong number of input arguments for rad_dke_yp');
0009 end
0010 %
0011 % Initialization of the coefficients
0012 %
0013 % Terms from Sp at l+1/2 (Notations: _imj -> (i,j+1+2);_ipj -> (i+1,j+1/2)
0014 %
0015 ZXXD_a.pp_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0016 ZXXD_a.pp_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0017 ZXXD_a.pm_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0018 ZXXD_a.pm_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0019 ZXXD_a.pr_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0020 ZXXD_a.pr_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0021 ZXXF_a.p_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0022 ZXXF_a.p_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0023 %
0024 % Terms from Sksi at l+1/2 (Notations: _ijp -> (i+1/2,j+1);_ijm -> (i+1/2,j)
0025 %
0026 ZXXD_a.mm_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0027 ZXXD_a.mm_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0028 ZXXD_a.mp_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0029 ZXXD_a.mp_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0030 ZXXD_a.mr_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0031 ZXXD_a.mr_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0032 ZXXF_a.m_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0033 ZXXF_a.m_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0034 %
0035 % Terms from Spsi at (i+1/2,j+1/2) (Notations: _lm -> (l);_lp -> (l+1);
0036 %
0037 ZXXD_a.rr_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0038 ZXXD_a.rr_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0039 ZXXD_a.rp_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0040 ZXXD_a.rp_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0041 ZXXD_a.rm_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0042 ZXXD_a.rm_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0043 ZXXF_a.r_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0044 ZXXF_a.r_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0045 %
0046 % Terms used by the old scheme only for momentum cross-derivatives (M.
0047 % Shoucri and I. Shkarofsky,Comp. Phys. Comm., 82 (1994) 287)
0048 %
0049 ZXXD_a.pp_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0050 ZXXD_a.pm_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0051 ZXXD_a.mp_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0052 ZXXF_a.p_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0053 %
0054 if ~isempty(transpfaste) && isstruct(transpfaste) && ~isempty(fieldnames(transpfaste)), 
0055     %
0056     [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0057     %
0058     npn = gridDKE.npn;
0059     nmhu = gridDKE.nmhu;
0060     nr = gridDKE.nr;
0061     xdpsin = gridDKE.xdpsin_dke;
0062     xdpsinm = gridDKE.xdpsinm_dke;
0063     xdpsinp = gridDKE.xdpsinp_dke;
0064     XXpn2 = gridDKE.XXpn2;
0065     Xmhu = gridDKE.Xmhu;
0066     Xmhu2 = gridDKE.Xmhu2;
0067     xmhubounce2 = gridDKE.xmhubounce2;
0068     pn = gridDKE.pn;
0069     ir_display = gridDKE.ir_display_out;
0070     mhu = gridDKE.mhu;
0071     %
0072     gamma = Zmomcoef.gamma;
0073     Xv = Zmomcoef.Xv;
0074     Xvpar = Zmomcoef.Xvpar;
0075     %
0076     ns_f = Zbouncecoef.ns_f;
0077     XXlambda = Zbouncecoef.XXlambda;
0078     blocksize_f_t = Zbouncecoef.blocksize_f_t;
0079     xqtilde = Zbouncecoef.xqtilde;
0080     %
0081     xx0 = equilDKE.xx0;
0082     xBp0 = equilDKE.xBp0;
0083     xB0 = equilDKE.xB0;
0084     ap = equilDKE.ap;
0085     Rp = equilDKE.Rp;
0086     psia_apRp = equilDKE.psia_apRp;
0087     xrho = equilDKE.xrho;
0088     %
0089     betath_ref = mksa.betath_ref;
0090     nhu_ref = mksa.nhu_ref;
0091     %
0092     bounce_mode = dkeparam.bounce_mode;
0093     %
0094     dp_cyl = display.dp_cyl;
0095     nlevel = display.nlevel;
0096     display_mode = display.display_mode;
0097     %
0098     time0 = clock;
0099     %
0100     [xTe_norm,xne_norm,xzni_norm,xzTi_norm,xnloss_norm,xbetath] = normcoef_dke_yp(mksa,equilDKE,gridDKE);   
0101     %
0102     XXDpsi = zeros(npn,nmhu,nr);%Defined on the distribution function grid (half grid) (psi coordinate)
0103     XXVpsi = zeros(npn,nmhu,nr);%Defined on the distribution function grid (half grid) (psi coordinate)
0104     %
0105     if transpfaste.norm_ref == 1,
0106         xvparmin_c = transpfaste.vparmin*ones(size(xrho));%Normalized to vth_ref
0107     elseif transpfaste.norm_ref == 0,
0108         xvparmin_c = transpfaste.vparmin*ones(size(xrho)).*xbetath/betath_ref;%Normalized to vth (local)
0109     end
0110     %
0111     xvparmin_r = real(xvparmin_c);
0112     xvparmin_i = imag(xvparmin_c);
0113     %
0114     % Various profile types
0115     %
0116     if isreal(transpfaste.pDr),
0117         xDr = real(transpfaste.Dr0)*(1 + transpfaste.pDr*xrho.^2);%flat or parabolic profile
0118     else
0119         xDr = (real(transpfaste.Dr0) - imag(transpfaste.Dr0))./(1+exp((xrho-real(transpfaste.pDr))/imag(transpfaste.pDr)))+imag(transpfaste.Dr0);%step profile
0120     end
0121     %
0122     if isreal(transpfaste.pVr),
0123         xVr = real(transpfaste.Vr0)*(1 + transpfaste.pVr*xrho.^2);%flat or parabolic profile
0124     else
0125         xVr = (real(transpfaste.Vr0) - imag(transpfaste.Vr0))./(1+exp((xrho-real(transpfaste.pVr))/imag(transpfaste.pVr)))+imag(transpfaste.Vr0);%step profile
0126     end
0127     %
0128     xmaxDpsi = zeros(1,nr);
0129     xmaxVpsi = zeros(1,nr);
0130     xEc_r_threshold = zeros(1,nr);
0131 %
0132     for ir = 1:nr,    
0133         ivpar_r_threshold = find(abs(Xvpar) >= xvparmin_r(ir),1);
0134         if isempty(ivpar_r_threshold),% xvparmin_r larger than last grid point
0135           ivpar_r_threshold = npn;
0136         end
0137         %
0138         xEc_r_threshold(ir) = mc2*(sqrt(1+pn(ivpar_r_threshold).*pn(ivpar_r_threshold)*betath_ref^2) - 1);%in keV
0139         %
0140         if real(transpfaste.Dr_model)==0,%No vpar dependence
0141             XDpsiv = xDr(ir)*(abs(Xvpar) >= xvparmin_r(ir))/(ap*ap*nhu_ref);%Normalized to the central collision frequency
0142         elseif real(transpfaste.Dr_model)==1,%Magnetic turbulence model: Dr(r)*vpar_ref/vth_ref
0143             XDpsiv = xDr(ir)*(abs(Xvpar) - xvparmin_r(ir)).*(abs(Xvpar) >= xvparmin_r(ir))/xvparmin_r(ir)/(ap*ap*nhu_ref);%Normalized to the central collision frequency
0144         elseif real(transpfaste.Dr_model)==-1,%Dr(r)*vth_ref/vpar_ref
0145             XDpsiv = xDr(ir)*(xvparmin_r(ir)*(abs(Xvpar) >= xvparmin_r(ir))./abs(Xvpar))/(ap*ap*nhu_ref);%Normalized to the central collision frequency
0146         elseif real(transpfaste.Dr_model)==2,
0147             if imag(transpfaste.vparmin) == 0,%No energy dependence
0148                 XDpsiv = xDr(ir)*(Xv >= xvparmin_r(ir))/(ap*ap*nhu_ref);%Normalized to the central collision frequency
0149             else
0150                 XDpsiv = xDr(ir)*(Xv >= (xvparmin_r(ir) - 2*xvparmin_i(ir))).*(Xv <= (xvparmin_r(ir) + 2*xvparmin_i(ir))).*exp(-(Xv - xvparmin_r(ir)).^2/xvparmin_i(ir).^2)/(ap*ap*nhu_ref);%Gaussian velocity dependence
0151             end
0152         end
0153         %
0154         XXDpsi(:,:,ir) = XDpsiv;%3-D radial diffusion matrix (radial profile dependence) (psi coordinate)
0155         xmaxDpsi(ir) = max(max(XXDpsi(:,:,ir)));%Maximum value of the radial diffusion coefficient (psi coordinate)
0156         %
0157         if real(transpfaste.Vr_model)==0,%No vpar dependence
0158             XVpsiv = xVr(ir)*(abs(Xvpar) >= xvparmin_r(ir))/(ap*nhu_ref);%Normalized to the central collision frequency (psi coordinate)
0159         elseif real(transpfaste.Vr_model)==1,%Magnetic turbulence model: Vr(r) = Vr0(r)*vpar_ref/vth_ref
0160             XVpsiv = xVr(ir)*(abs(Xvpar) - xvparmin_r(ir)).*(abs(Xvpar) >= xvparmin_r(ir))/xvparmin_r(ir)/(ap*nhu_ref);%Normalized to the central collision frequency (psi coordinate)
0161         elseif real(transpfaste.Vr_model)==-1,%Vr(r) = Vr0(r)*vth_ref/vpar_ref
0162             XVpsiv = xVr(ir)*(xvparmin_r(ir)*(abs(Xvpar) >= xvparmin_r(ir))./abs(Xvpar))/(ap*nhu_ref);%Normalized to the central collision frequency (psi coordinate)
0163         elseif real(transpfaste.Vr_model)==2,
0164             if imag(transpfaste.vparmin) == 0,%No energy dependence
0165                 XVpsiv = xVr(ir)*(Xv >= xvparmin_r(ir))/(ap*nhu_ref);%Normalized to the central collision frequency (psi coordinate)
0166             else
0167                 XVpsiv = xVr(ir)*exp(-(Xv - xvparmin_r(ir)).^2/xvparmin_i(ir).^2)/(ap*nhu_ref);
0168             end
0169         end
0170         %
0171         XXVpsi(:,:,ir) = XVpsiv;%3-D radial pinch matrix (radial profile dependence) (psi coordinate)
0172         xmaxVpsi(ir) = max(max(XXVpsi(:,:,ir)));%Maximum value of the radial convection coefficient (psi coordinate)
0173         %
0174     end
0175 %
0176     if display_mode == 2,
0177         Fig_radialdiff = figure('Name','Fast electron radial diffusion');
0178         %
0179         figure(Fig_radialdiff),subplot(121);zoom on;
0180         [XDpsiv_cyl,ppar_cyl,dppar_cyl,pperp_cyl] = s2c_dke_yp(XXDpsi(:,:,ir_display),pn,mhu,dp_cyl);
0181         %
0182         Xppar_cyl = ones(length(pperp_cyl),1)*ppar_cyl(:)';
0183         Xpperp_cyl = pperp_cyl(:)*ones(1,length(ppar_cyl));
0184         Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0185         %
0186         XDpsiv_cyl = XDpsiv_cyl.*(XDpsiv_cyl>0);%Remove negative values
0187         XDpsiv_cyl = XDpsiv_cyl.*(Xp_cyl<max(ppar_cyl));%Remove values above max(ppar)
0188         %
0189         contourf(ppar_cyl,pperp_cyl,XDpsiv_cyl,nlevel);
0190         if xmaxDpsi(ir_display)>0,colorbar;end
0191         axis('equal');
0192         faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0193         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,['rho = ',num2str(xrho(ir_display))]);
0194         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*10,['v//-threshold(norm) = ',num2str(xvparmin_r(ir_display))]);
0195         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*15,['Ec-threshold = ',num2str(xEc_r_threshold(ir_display)),' (keV)']);
0196         xlabel('ppar/pth-ref');ylabel('pperp/pth-ref');title('Radial diffusion coefficient (nhu-ref*psia_apRp*psia_apRp)');
0197         %
0198         subplot(122),plot(pn./gamma,XXDpsi(:,nmhu,ir_display),'b');
0199         faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0200         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,'mhu = 1');
0201         xlabel('v/vth-ref');ylabel('Dr (nhu-ref*psia_apRp*psia_apRp)');title('Radial diffusion coefficient ');drawnow
0202         %
0203         Fig_radialpinch = figure('Name','Fast electron radial pinch');
0204         figure(Fig_radialpinch),subplot(121);zoom on,
0205         [XVpsiv_cyl,ppar_cyl,dppar_cyl,pperp_cyl] = s2c_dke_yp(XXVpsi(:,:,ir_display),pn,mhu,dp_cyl);
0206         %
0207         Xppar_cyl = ones(length(pperp_cyl),1)*ppar_cyl(:)';
0208         Xpperp_cyl = pperp_cyl(:)*ones(1,length(ppar_cyl));
0209         Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0210         %
0211         XVpsiv_cyl = XVpsiv_cyl.*(XVpsiv_cyl>0);%Remove negative values
0212         XVpsiv_cyl = XVpsiv_cyl.*(Xp_cyl<max(ppar_cyl));%Remove values above max(ppar)
0213         %
0214         contourf(ppar_cyl,pperp_cyl,XVpsiv_cyl,nlevel);
0215         if xmaxVpsi(ir_display)>0,colorbar;end
0216         axis('equal');
0217         faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0218         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,['rho = ',num2str(xrho(ir_display))]);
0219         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*10,['v//-threshold(norm) = ',num2str(xvparmin_r(ir_display))]);
0220         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*15,['Ec-threshold = ',num2str(xEc_r_threshold(ir_display)),' (keV)']);
0221         xlabel('ppar/pth-ref');ylabel('pperp/pth-ref');title('Radial convection coefficient (nhu-ref*psia_apRp)');
0222         %
0223         subplot(122),plot(pn./gamma,XXVpsi(:,nmhu,ir_display),'r');
0224         faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0225         text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,'mhu = 1');
0226         xlabel('v/vth-ref');ylabel('Vr (nhu-ref*psia_apRp)');title('Radial convection coefficient ');drawnow
0227         %
0228         Fig_radialtranspvthreshold = figure('Name','Fast electron radial transport velocity threshold');
0229         figure(Fig_radialtranspvthreshold);
0230         subplot(121),plot(xrho,xvparmin_r,'g-x');grid on;zoom on;
0231         axis([0,1,0,1.1*max(xvparmin_r)]);
0232         xlabel('rho');
0233         ylabel('v// threshold (vth-ref)');
0234         title('Radial transport');
0235         subplot(122),plot(xrho,xEc_r_threshold,'c-x');grid on;zoom on;
0236         axis([0,1,0,1.1*max(xEc_r_threshold)]);
0237         xlabel('xrho');
0238         ylabel('Kinetic energy threshold (keV)');
0239         title('Radial transport');drawnow
0240         %
0241         Fig_fastetransp = figure('Name','Fast electron radial transport'); 
0242         figure(Fig_fastetransp);
0243         subplot(121),plot(xrho,xmaxDpsi,'b-o');grid on;zoom on;
0244         if sum(xmaxDpsi) ~= 0,axis([0,1,0,1.1*max(xmaxDpsi)]);end
0245         xlabel('rho');
0246         ylabel('Maximum value (nhu-ref*psia_apRp*psia_apRp)');
0247         title('Radial diffusion');
0248         subplot(122),plot(xrho,xmaxVpsi,'r-o');grid on;zoom on;
0249         if sum(xmaxVpsi) ~= 0,axis([0,1,0,1.1*max(xmaxVpsi)]);end
0250         xlabel('rho');
0251         ylabel('Maximum value (nhu-ref*psia_apRp)');
0252         title('Radial convection');drawnow
0253     end
0254     %
0255     xdpsinmm = [0,xdpsin(1:nr-1)];
0256     xdpsinpp = [xdpsin(2:nr),0];
0257     %
0258     XXdpsin = reshape(repmat(xdpsin,npn*nmhu,1),[npn,nmhu,nr]);   
0259     XXdpsinm = reshape(repmat(xdpsinm,npn*nmhu,1),[npn,nmhu,nr]); 
0260     XXdpsinp = reshape(repmat(xdpsinp,npn*nmhu,1),[npn,nmhu,nr]); 
0261     XXdpsinmm = reshape(repmat(xdpsinmm,npn*nmhu,1),[npn,nmhu,nr]);   
0262     XXdpsinpp = reshape(repmat(xdpsinpp,npn*nmhu,1),[npn,nmhu,nr]);
0263     %
0264     XXdeltalp = XXdpsinpp./(XXdpsinpp + XXdpsin);%Grid interpolation coefficients
0265     XXdeltalm = XXdpsin./(XXdpsinmm + XXdpsin);%Grid interpolation coefficients
0266     %
0267     XXDpsimm = zeros(npn,nmhu,nr);
0268     XXDpsimm(:,:,2:nr) = XXDpsi(:,:,1:nr-1);
0269     XXDpsipp = zeros(npn,nmhu,nr);
0270     XXDpsipp(:,:,1:nr-1) = XXDpsi(:,:,2:nr);
0271     %
0272     XXVpsimm = zeros(npn,nmhu,nr);
0273     XXVpsimm(:,:,2:nr) = XXVpsi(:,:,1:nr-1);
0274     XXVpsipp = zeros(npn,nmhu,nr);
0275     XXVpsipp(:,:,1:nr-1) = XXVpsi(:,:,2:nr);
0276     %
0277     % Terms from Spsi at (i+1/2,j+1/2) (Notations: _lm -> (l);_lp -> (l+1);
0278     %
0279     ZXXD_a.rr_lm = (1 - XXdeltalm).*XXDpsi + XXdeltalm.*XXDpsimm;%Evaluation of the radial diffusion rate on the flux grid (linear interpolation)
0280     ZXXD_a.rr_lp = (1 - XXdeltalp).*XXDpsipp + XXdeltalp.*XXDpsi;%Evaluation of the radial diffusion rate on the flux grid (linear interpolation)
0281     ZXXF_a.r_lm = (1 - XXdeltalm).*XXVpsi + XXdeltalm.*XXVpsimm;%Evaluation of the radial convection on the flux grid (linear interpolation)
0282     ZXXF_a.r_lp = (1 - XXdeltalp).*XXVpsipp + XXdeltalp.*XXVpsi;%Evaluation of the radial convection on the flux grid (linear interpolation)
0283 end
0284 
0285
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.