finit_dke_yp

   Calculate initial distribution function for the 3D electron relativistic drift kinetic solver.

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

0001 function [XXfinit] = finit_dke_yp(dkeparam,display,equilDKE,gridDKE,mksa,Zmomcoef,waveparam,XXf0_in)
0002 %
0003 %   Calculate initial distribution function for the 3D electron relativistic drift kinetic solver.
0004 %
0005 %   by Y. Peysson (CEA-DRFC) (yves.peysson@cea.fr) and J. Decker (CEA-DRFC) (joan.decker@cea.fr)
0006 %
0007 npn = gridDKE.npn;
0008 nmhu = gridDKE.nmhu;
0009 nr = gridDKE.nr;
0010 ir_display = gridDKE.ir_display_out;
0011 pn = gridDKE.pn;
0012 pn2 = gridDKE.pn2;
0013 dpn = gridDKE.dpn;
0014 mhu = gridDKE.mhu;
0015 dmhu = gridDKE.dmhu;
0016 Xpn = gridDKE.Xpn;
0017 Xmhu = gridDKE.Xmhu;
0018 zZi = equilDKE.zZi;
0019 %
0020 xrho = equilDKE.xrho;
0021 %
0022 gamma = Zmomcoef.gamma;
0023 Xgamma = Zmomcoef.Xgamma;
0024 v = Zmomcoef.v;
0025 %
0026 if isfield(display,'display_init'),
0027     display_mode = display.display_init;
0028 else
0029     display_mode = 0;
0030 end
0031 %
0032 dp_cyl = display.dp_cyl;
0033 nlevel = display.nlevel;
0034 %
0035 dke_mode = dkeparam.dke_mode;
0036 %
0037 [xTe_norm,xne_norm,xzni_norm,xzTi_norm,xnloss_norm,xbetath,xlnc_e,xnhu,xrnhuth] = normcoef_dke_yp(mksa,equilDKE,gridDKE);
0038 %
0039 XXfinit = zeros(npn,nmhu,nr);%Defined on the distribution function grid (half grid)
0040 %
0041 if dkeparam.initguess_f == 1 && exist('waveparam'),
0042     if isreal(waveparam.bdNpar_rf),%General RF power spectrum
0043          dkeparam.initguess_f = 0;
0044     else,%dNpar is made imaginary to identify square power spectrum type (for the Lower Hybrid wave comparison with theoretical models only)
0045          xvparmin_lh = 1./(waveparam.bNpar_rf + imag(waveparam.bdNpar_rf)/2)/mksa.betath_ref;
0046          xvparmax_lh = 1./(waveparam.bNpar_rf - imag(waveparam.bdNpar_rf)/2)/mksa.betath_ref;
0047     end
0048 else
0049     dkeparam.initguess_f = 0;%Maxwellian distribition function as initial guess if no RF wave in the plasma
0050 end
0051 %
0052 for ir = 1:nr,    
0053     %
0054     fM = maxwellian_dke_yp(gamma(:),pn2(:),dpn(:),xne_norm,xTe_norm,ir);%Maxwellian distribution normalized to normalized density ne_norm (1st order gamma corrections)
0055     XfM = fM*ones(1,nmhu);
0056     %
0057     if dkeparam.initguess_f == 0 && isempty(XXf0_in),
0058         Xf0 = XfM;
0059         title_initguess = 'Relativistic Maxwellian';        
0060     elseif dkeparam.initguess_f == 1 && isempty(XXf0_in),
0061 %
0062 %   Initial guess for the distribution function from ref. V.Fuchs et al.,Phys. Fluids 28 (1985) 3619.
0063 %   (Only for the old method for calculating the quasilinear diffusion coefficient of the LH operator)
0064 %
0065         alpha = 2 - 2/(2+zZi(1));%Dominant ion contribution (it is a guess only)
0066         beta = 2 - alpha;
0067         Tp1nr = (xvparmax_lh^(alpha+1) - xvparmin_lh^(alpha+1))*(alpha-1)/(alpha+1) - (alpha-1)*(xvparmax_lh - xvparmin_lh)*xvparmin_lh^alpha + xvparmax_lh^(alpha-1) - xvparmin_lh^(alpha-1);
0068         Tp2nr = 2*(alpha-1)*(xvparmax_lh - xvparmin_lh)*xvparmin_lh^(alpha-2) - (xvparmax_lh^(alpha-1) - xvparmin_lh^(alpha-1));
0069         Tpnr = Tp1nr/Tp2nr;
0070         Tpr = (1 + 8*mksa.betath_ref/5)*Tpnr;%Relativistic correction
0071         Xpperp2 = Xpn.*Xpn.*(1 - Xmhu.*Xmhu);
0072         Xppar = Xpn.*Xmhu;
0073         Xppar2 = Xppar.*Xppar;
0074         Xpparmin = xvparmin_lh*sqrt((1+Xpperp2*mksa.betath_ref^2)/(1-mksa.betath_ref^2*xvparmin_lh^2));
0075         Xpparmax = xvparmax_lh*sqrt((1+Xpperp2*mksa.betath_ref^2)/(1-mksa.betath_ref^2*xvparmax_lh^2));
0076         Xmaskplat = (Xppar >= Xpparmin).*(Xppar <= Xpparmax);
0077         Xfsth0 = Xmaskplat.*exp(-Xpperp2/2/Tpr)/2/pi/Tpr;
0078         Xmaskback = (Xppar <= Xpparmin);
0079         Xfsth1 = Xmaskback.*exp(-Xpn.*Xpn/2/Tpr).*exp(Xpparmin.*Xpparmin/2/Tpr)/2/pi/Tpr;
0080         Xmaskup = (Xppar > Xpparmax);
0081         Xfsth2 = Xmaskup.*exp(-(Xpn.*Xpn.*(Xppar./Xpparmax).^(2/zZi(1)) - Xpparmax.*Xpparmax)/Tpr/2)/2/pi/Tpr;%See thesis of R.Arslanbekov
0082         Xfsth = Xfsth0+Xfsth1+Xfsth2;
0083         fM = xne_norm(ir)*exp(-pn.*pn./(1+gamma)/xTe_norm(ir))/2/pi/xTe_norm(ir)/sqrt(2*pi*xTe_norm(ir));
0084         XfM = fM'*ones(1,nmhu);
0085         normfM = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,XfM,2),1);
0086         normfsth = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xfsth,2),1);
0087         XfM = XfM/normfM;
0088         Xfsth = Xfsth/normfsth;
0089         lambda = 1/(1 + Xfsth(min(find(Xmaskplat(:,nmhu-1)==1)),nmhu-1)/XfM(min(find(Xmaskplat(:,nmhu-1)==1)),nmhu-1));
0090         Xf0 = XfM + Xfsth*lambda;   
0091         normf0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xf0,2),1);
0092         Xf0 = xne_norm(ir)*Xf0/normf0;
0093         title_initguess = 'V. Fuchs et al., Phys. Fluids, 28 (1985) 3619';
0094     else
0095         if isstruct(XXf0_in),
0096             if ~comp_struct_jd(mksa,XXf0_in.mksa,0),
0097                 error('initial distribution cannot be used with a different normalization. Use dkeparam.Te_ref and dkeparam.ne_ref, or use loadstructs')
0098             end
0099             if size(XXf0_in.XXf0_leg,2) == nmhu,
0100                 Xf0 = XXf0_in.XXf0_leg(:,:,ir);
0101                 title_initguess = 'Initial distribution from input arguments';
0102             else
0103                 Xf0 = interpleg_yp(XXf0_in.XXf0_leg(:,:,ir),mhu);%interpolate using Legendre polynomials
0104                 title_initguess = 'Initial distribution interpolated from input arguments';
0105             end
0106             %
0107             if isfield(dkeparam,'ne_ref'),
0108                 %
0109                 % density variations between time steps are accounted for by adding or substracting a maxwellian component
0110                 %
0111                 normf0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xf0,2),1);
0112                 Xf0 = Xf0 + (1 - normf0/xne_norm(ir))*XfM;
0113                 %
0114                 % avoiding negative values for the distribution
0115                 %
0116                 Xf0(Xf0 < 0) = 0;
0117                 normf0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xf0,2),1);
0118                 Xf0 = Xf0*xne_norm(ir)/normf0;
0119             end
0120             %
0121         else
0122             if size(XXf0_in,2) == nmhu,
0123                 Xf0 = XXf0_in(:,:,ir);
0124                 title_initguess = 'Initial distribution from input arguments';
0125             else
0126                 Xf0 = interpleg_yp(XXf0_in(:,:,ir),mhu);%interpolate using Legendre polynomials
0127                 title_initguess = 'Initial distribution interpolated from input arguments';
0128             end
0129         end
0130     end
0131     %
0132     XXfinit(:,:,ir) = Xf0;
0133 end
0134 %
0135 % *************************** Display initial distribution function ***************************
0136 %
0137 if display_mode ~= 0,
0138     Fig_InitDistr = figure('Name','Initial distribution');    
0139     Fig_InitDistrp0 = figure('Name','Initial distribution at pperp=0');
0140     Fig_ElecVeloc = figure('Name','Electron velocity');
0141 end
0142 %
0143 for ir = 1:nr,
0144     if display_mode ~= 0,
0145         if dke_mode == 1,
0146 %
0147             view = 0;
0148 %
0149             if ir == ir_display-1,
0150                 figure(Fig_InitDistr),subplot(131);zoom on
0151                 view = 1;
0152             elseif ir == ir_display,
0153                 figure(Fig_InitDistr),subplot(132);zoom on
0154                 view = 1;
0155             elseif ir == ir_display+1,
0156                 figure(Fig_InitDistr),subplot(133);zoom on
0157                 view = 1;
0158             end    
0159 %
0160             if view == 1,        
0161                 normf0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xf0,2),1);
0162                 [logXf0_cyl,ppar_cyl,dppar_cyl,pperp_cyl,dpperp_cyl] = s2c_dke_yp(log(Xf0),pn,mhu,dp_cyl);
0163                 Xf0_cyl = exp(logXf0_cyl);%For accurate representation in figures
0164 %
0165                 Xppar_cyl = ones(length(pperp_cyl),1)*ppar_cyl(:)';
0166                 Xpperp_cyl = pperp_cyl(:)*ones(1,length(ppar_cyl));
0167                 Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0168                 Xmhu_cyl = Xppar_cyl./sqrt(Xppar_cyl.^2+Xpperp_cyl.^2);
0169                 Xgamma_cyl = sqrt(1 + Xp_cyl.*Xp_cyl*mksa.betath_ref*mksa.betath_ref);
0170                 Xvpar_cyl = Xppar_cyl./Xgamma_cyl;
0171 %
0172                 Xf0_cyl = Xf0_cyl.*(Xf0_cyl>0);%Remove negative values
0173                 Xf0_cyl = Xf0_cyl.*(Xp_cyl<max(ppar_cyl));%Remove values above max(ppar)
0174                 Xf0_cyl(isnan(Xf0_cyl)) = 0;%Remove NaN values
0175 %
0176                 curr0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*(pn(:)./gamma(:)).*integral_dke_jd(dmhu,Xmhu.*Xf0,2),1);
0177                 curr0_cyl = 2*pi*integral_dke_jd(dpperp_cyl(:),pperp_cyl(:).*integral_dke_jd(dppar_cyl,Xvpar_cyl.*Xf0_cyl,2),1);
0178 %
0179                 C = contour(ppar_cyl,pperp_cyl,log(Xf0_cyl),[-20:20/nlevel:0]);axis('equal');
0180                 
0181                 
0182                 
0183                 
0184                 
0185                 
0186                 faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0187                 text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*5,['rho = ',num2str(xrho(ir))]);
0188                 text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*10,['Te(norm) = ',num2str(xTe_norm(ir))]);
0189                 text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*15,['ne(norm) = ',num2str(xne_norm(ir))]);
0190                 xlabel('ppar/pth-ref');ylabel('pperp/pth-ref');title(title_initguess);
0191                 drawnow;
0192 %
0193             end
0194 %
0195             view = 0;
0196             if ir == ir_display-1,
0197                 figure(Fig_InitDistrp0),subplot(231);zoom on
0198                 view = 1;
0199             elseif ir == ir_display,
0200                 figure(Fig_InitDistrp0),subplot(232);zoom on
0201                 view = 1;
0202             elseif ir == ir_display+1,
0203                 figure(Fig_InitDistrp0),subplot(233);zoom on
0204                 view = 1;
0205             end
0206 %
0207             if view == 1,
0208                 semilogy(-v,XfM(:,1),'g--',v,XfM(:,nmhu),'g--',-v,Xf0(:,1),'c-',v,Xf0(:,nmhu),'c-');%Forward and backward directions only
0209                 axis([-max(v),max(v),1e-20,10]);grid on;zoom on
0210                 xlabel('vpar/vth-ref');ylabel('f(pperp = 0)');title(['Norm.(rho=',num2str(xrho(ir)),') = ',num2str(normf0)]);
0211                 drawnow;
0212             end
0213 %
0214             view = 0;
0215             if ir == ir_display-1,
0216                 figure(Fig_ElecVeloc),subplot(234);zoom on
0217                 view = 1;
0218             elseif ir == ir_display,
0219                 figure(Fig_ElecVeloc),subplot(235);zoom on
0220                 view = 1;
0221             elseif ir == ir_display+1,
0222                 figure(Fig_ElecVeloc),subplot(236);zoom on
0223                 view = 1;
0224             end  
0225 %
0226             if view == 1,
0227                 plot(pn,pn,'r-.',pn,v,'b:');grid on;zoom on
0228                 xlabel('p/pth-ref');ylabel('v/vth-ref');title(['j = ',num2str(curr0)]);
0229                 drawnow;
0230             end
0231         else
0232             view = 0;
0233             if ir == ir_display,
0234                 figure(Fig_InitDistr),zoom on
0235                 view = 1;
0236             end
0237 %
0238             if view == 1,
0239                 normf0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*integral_dke_jd(dmhu,Xf0,2),1);
0240                 [logXf0_cyl,ppar_cyl,dppar_cyl,pperp_cyl,dpperp_cyl] = s2c_dke_yp(log(Xf0),pn,mhu,dp_cyl);
0241                 Xf0_cyl = exp(logXf0_cyl);%For accurate representation in figures
0242 %
0243                 Xppar_cyl = ones(length(pperp_cyl),1)*ppar_cyl(:)';
0244                 Xpperp_cyl = pperp_cyl(:)*ones(1,length(ppar_cyl));
0245                 Xp_cyl = sqrt(Xppar_cyl.*Xppar_cyl + Xpperp_cyl.*Xpperp_cyl);
0246                 Xmhu_cyl = Xppar_cyl./sqrt(Xppar_cyl.^2 + Xpperp_cyl.^2);
0247                 Xgamma_cyl = sqrt(1 + Xp_cyl.*Xp_cyl*mksa.betath_ref*mksa.betath_ref);
0248                 Xvpar_cyl = Xppar_cyl./Xgamma_cyl;
0249 %
0250                 Xf0_cyl = Xf0_cyl.*(Xf0_cyl>0);%Remove negative values
0251                 Xf0_cyl = Xf0_cyl.*(Xp_cyl<max(ppar_cyl));%Remove values above max(ppar)
0252 %
0253                 curr0 = 2*pi*integral_dke_jd(dpn(:),pn2(:).*(pn(:)./gamma(:)).*integral_dke_jd(dmhu,Xmhu.*Xf0,2),1);
0254                 curr0_cyl = 2*pi*integral_dke_jd(dpperp_cyl(:),pperp_cyl(:).*integral_dke_jd(dppar_cyl,Xvpar_cyl.*Xf0_cyl,2),1);
0255 %
0256                 graph2D_jd(ppar_cyl,pperp_cyl,Xf0_cyl','p_{||}/p_{Te}','p_{\perp}/p_{Te}',title_initguess,[-dkeparam.pnmax_S,dkeparam.pnmax_S],[0,dkeparam.pnmax_S],0,mksa.betath_ref,[0:1:dkeparam.pnmax_S],0,'-','r',2,20,0.9,max(max(Xf0_cyl)));
0257                 axis('equal');axis([-dkeparam.pnmax_S,dkeparam.pnmax_S,0,dkeparam.pnmax_S]);drawnow;
0258                 faxis = axis;dfaxis_x = (faxis(2) - faxis(1))/100;dfaxis_y = (faxis(4) - faxis(3))/100; 
0259                 text(faxis(1) + dfaxis_x*5,faxis(4) - dfaxis_y*6,['\rho = ',num2str(xrho(ir))]);
0260 %
0261             end
0262 %
0263             view = 0;
0264             if ir == ir_display,
0265                 figure(Fig_InitDistrp0),zoom on;drawnow;
0266                 view = 1;
0267             end
0268 %
0269             if view == 1,         
0270                 graph1D_jd(-v,XfM(:,1),0,1,'','','',NaN,[-max(v),max(v)],[1e-20,10],'-','none','b',1,20,gca,0.9,0.7,0.7);
0271                 graph1D_jd(v,XfM(:,nmhu),0,1,'','','',NaN,[-max(v),max(v)],[1e-20,10],'-','none','b',1,20,gca,0.9,0.7,0.7);
0272                 graph1D_jd(-v,Xf0(:,1),0,1,'','','',NaN,[-max(v),max(v)],[1e-20,10],'-','none','r',2,20,gca,0.9,0.7,0.7);
0273                 graph1D_jd(v,Xf0(:,nmhu),0,1,'v_{||}/v_{Te}','f(p_{\perp}=0)',['Norm.=',num2str(normf0),', \beta_{Te}=',num2str(mksa.betath_ref),', J_{LH}=',num2str(curr0)],NaN,[-max(v),max(v)],[1e-20,10],'-','none','r',2,20,gca,0.9,0.7,0.7);
0274                 grid on;zoom on;drawnow;                
0275             end
0276 %
0277             view = 0;
0278             if ir == ir_display,
0279                 figure(Fig_ElecVeloc);zoom on;drawnow;
0280                 view = 1;
0281             end    
0282 %
0283             if view == 1,
0284                 graph1D_jd(pn,pn,0,0,'','','',NaN,[0,max(pn)],[0,max(v)],'-','none','b',1,20,gca,0.9,0.7,0.7);
0285                 graph1D_jd(pn,v,0,0,'p/p_{Te}','v/v_{Te}','',NaN,[0,max(pn)],[0,max(v)],'-','none','r',2,20,gca,0.9,0.7,0.7);
0286                 grid on;zoom on;drawnow;                
0287             end
0288 %
0289         end
0290     end
0291 end