mripplecoefbuilder_dke_yp

   Calculate magnetic ripple coefficients 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 [Zlossmripple] = mripplecoefbuilder_dke_yp(dkeparam,display,equil,equilDKE,gridDKE,mksa,ripple)
0002 %
0003 %   Calculate magnetic ripple coefficients 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 %
0008 npn = gridDKE.npn;
0009 nmhu = gridDKE.nmhu;
0010 nr = gridDKE.nr;
0011 xrho = gridDKE.xrho; 
0012 %
0013 display_mode = display.display_mode;
0014 %
0015 [xTe_norm,xne_norm,xzni_norm,xzTi_norm,xnloss_norm,xbetath,xlnc_e,xnhu,xrnhuth] = normcoef_dke_yp(mksa,equilDKE,gridDKE);
0016 %
0017 XXlossripple = zeros(npn,nmhu,nr);%Defined on the distribution function grid (half grid)
0018 xdeltaripple = zeros(size(xrho));
0019 xdeltaripple = zeros(size(xrho));
0020 xNq = zeros(size(xrho));
0021 Xnhuloss = zeros(npn,nmhu);
0022 xthetabanana = zeros(size(xrho));
0023 xthetaloss = zeros(size(xrho));
0024 xpndetrap = zeros(size(xrho));
0025 xEseuil = zeros(size(xrho));
0026 %
0027 if ~isempty(ripple) && (isstruct(ripple) && isfield(ripple,'mode') && ripple.mode == 1),
0028     etime_rip = 0;
0029     time0 = clock;
0030     %
0031     bounce_mode = dkeparam.bounce_mode;
0032     %
0033     xne = equilDKE.xne;
0034     xBT0 = equilDKE.xBT0;
0035     xZeff = equilDKE.xZeff;
0036     ap = equilDKE.ap;
0037     Rp = equilDKE.Rp;
0038     %
0039     Xpn = gridDKE.Xpn;
0040     Xmhu = gridDKE.Xmhu;
0041     %
0042     tokname = equil.id;%Tokamak name
0043     %
0044     xnhuloss = ones(size(xrho))*1000;%Normalized to the reference frequency(forced to very large loss frequency, otherwise no physical meaning)
0045     %
0046     for ir = 1:nr,
0047         if isnan(xnhuloss(ir)) == 0 && bounce_mode == 1 && ~(isempty(findstr(tokname,'TS')) && isempty(findstr(tokname,'WEST'))) == 1,%WARNING: only valid for the TORE SUPRA tokamak with large magnetic ripple
0048             if ir == 1, 
0049                 info_dke_yp(2,'WARNING: Ripple losses calculations only valid for Tore Supra and WEST tokamak, see main_dke_yp.m for BeLi, Ip and Ncoils values');
0050                 disp('-------------------------------------------------------------------------------------------------------------------');
0051                 Ip = 0.61;%MA (For magnetic equilibrium in the magnetic ripple routine)
0052                 BeLi = 1.0;%For magnetic equilibrium in the magnetic ripple routine
0053             end
0054 %            Ncoils = 1;%Ncoils is here forced to 1, as prescribed in Ref. M. Ju et al., Phys. Plasmas, 9 (2002) 493. Upper bound of ripple losses is then calculated
0055 %            xNq(ir) = Ncoils*xq(ir);%Exact N may be used but large values may lead to numerical instabilities due to round-off errors
0056             [xthetabanana(ir),xthetaloss(ir),pndetrap,xEseuil(ir),xdeltaripple(ir)] = detrapmj3(xZeff(ir),ap,xBT0(ir),Rp,BeLi,Ip,xne(ir),xrho(ir));%WARNING: Only valid for the TORE SUPRA tokamak
0057             xpndetrap(ir) = pndetrap/mksa.betath_ref;
0058 %
0059             xmhusupertrap(ir) = sin(pi*xthetaloss(ir)/180);
0060             Xmaskrip = (Xpn > xpndetrap(ir)).*(abs(Xmhu) < xmhusupertrap(ir));
0061             Xnhuloss = Xmaskrip*xnhuloss(ir);
0062         else
0063             if isnan(xnhuloss(ir)) == 0 & bounce_mode == 1 & (isempty(findstr(tokname,'TS')) || isempty(findstr(tokname,'WEST'))) == 1
0064                 if ir == 1, 
0065                     info_dke_yp(2,'WARNING: Ripple losses calculations only valid for Tore Supra or WEST tokamak. Wrong tokamak name');
0066                     disp('-------------------------------------------------------------------------------------------------------------------');
0067                 end
0068             elseif isnan(xnhuloss(ir)) == 0 & bounce_mode == 0 & (isempty(findstr(tokname,'TS')) || isempty(findstr(tokname,'WEST'))) == 0
0069                 if ir == 1, 
0070                     info_dke_yp(2,'WARNING: Ripple losses calculations only valid for bounce averaged calculations');
0071                     disp('-------------------------------------------------------------------------------------------------------------------');
0072                 end
0073             end    
0074         end
0075         XXlossripple(:,:,ir) = Xnhuloss;%Magnetic ripple loss term is expressed in a Krook term
0076     end
0077     %
0078     etime_rip = etime_rip + etime(clock,time0);
0079     if display_mode >= 1,info_dke_yp(2,['Magnetic ripple coefficients calculations done in ',num2str(etime_rip),' (s)']);end    
0080 end
0081 %
0082 %if isempty(XXsource),XXsource = zeros(npn,nmhu,nr);end
0083 %XXsource = XXsource/mksa.ne_ref/mksa.nhu_ref;%Source or sink terms
0084 %
0085 Zlossmripple.XXlossripple = XXlossripple;
0086 Zlossmripple.xdeltaripple = xdeltaripple;
0087 Zlossmripple.xNq = xNq;
0088 Zlossmripple.Xnhuloss = Xnhuloss;
0089 Zlossmripple.xthetabanana = xthetabanana;
0090 Zlossmripple.xthetaloss = xthetaloss;
0091 Zlossmripple.xpndetrap = xpndetrap;
0092 Zlossmripple.xpnth = xbetath./xbetath(1);
0093 Zlossmripple.xEseuil = xEseuil;
0094 %
0095