loop_run_lukert

 Ray-tracing + quasilinear kinetic calculations using LUKE. Tool for implementation in other calculation packages.

 by Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr> and J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr>

0001 function [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = loop_run_lukert(simul,waves,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct)
0002 %
0003 % Ray-tracing + quasilinear kinetic calculations using LUKE. Tool for implementation in other calculation packages.
0004 %
0005 % by Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr> and J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr>
0006 %
0007 path_simul = '';
0008 %
0009 nw = length(waves);
0010 nb = zeros(1,nw);
0011 for iw = 1:nw,
0012     if isfield(waves{iw},'rays'),
0013         nb(iw) = length(waves{iw}.rays);
0014     else
0015         nb(iw) = NaN;
0016     end
0017 end
0018 %
0019 % Determination of the radial grid optmized according to linear damping
0020 %
0021 if ~isfield(dkeparam,'w_mask'),
0022     w_mask = NaN;
0023 else
0024     w_mask = dkeparam.w_mask;%list of waves to account for in the radial grid determination
0025 end
0026 %
0027 rho_S_init = dkeparam.rho_S;
0028 %
0029 if isfield(dkeparam,'rho0') && length(rho_S_init) == 1,
0030     rho_S_init = rho_S_init + 1i*dkeparam.rho0;
0031 end
0032 %
0033 if isfield(dkeparam,'nk_abs'),
0034     nk_abs = dkeparam.nk_abs;
0035 else
0036     nk_abs = 1;
0037 end
0038 %
0039 [dkeparam.rho_S,dkeparam.pnmax_S(1)] = optgrid_dke_jd(waves,rho_S_init,dkeparam.pnmax_S(1),equil,w_mask,dkedisplay.display_mode);
0040 %
0041 % 3-D relativistic bounce averaged quasilinear kinetic calculations (LUKE)
0042 %
0043 if dkeparam.rt_mode == 0 || dkeparam.rt_mode == 1,
0044 %
0045 % Case with sdNpar given by the waves structure (usualy constant along ray
0046 % propagation) if the waves structure is filled, or case without RF power
0047 % if waves is empty
0048 %
0049     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0050 %
0051 % Case with sdNpar proportional to sNpar (in waves structure)
0052 %
0053 elseif dkeparam.rt_mode == 2,
0054     %
0055     [qe,me] = pc_dke_yp;
0056     %
0057     for iw = 1:nw,
0058         for ib = 1:nb(iw),
0059             %
0060             dNpar0 = waves{iw}.rays{ib}.sdNpar(1);
0061             Npar0 = real(waves{iw}.rays{ib}.sNpar(1));
0062             %
0063             npresmin = 5;%minimum number of resonant points for numerical stability
0064             vparres = 4;%resonant velocity
0065             %
0066             dNparmin = abs(Npar0)*npresmin*dkeparam.pnmax_S/dkeparam.np_S/vparres;
0067             %
0068             omega_ce = qe*equil.ptBPHI(1,1)/me;
0069             %
0070             if waves{iw}.omega_rf < omega_ce/5 && dNpar0 < dNparmin,
0071                 info_dke_yp(2,['WARNING: wave #',num2str(iw),' ray #',num2str(ib),': initial dNpar increased to ',num2str(dNparmin),' for numerical stability']);
0072                 dNpar0 = dNparmin;
0073             end
0074             %
0075             waves{iw}.rays{ib}.sdNpar = real(waves{iw}.rays{ib}.sNpar)*dNpar0/Npar0;
0076         end
0077     end
0078     %
0079     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0080 %
0081 % *************************************************************************
0082 %
0083 % Case with sdNpar eventually proportional to ss (in waves structure)
0084 %
0085 elseif dkeparam.rt_mode == 4,
0086     %
0087     iw_list = [];%list of LH waves
0088     a_sdNpar = zeros(1,nw);
0089     %
0090     for iw = 1:nw,
0091         if strcmp(waves{iw}.rayinit.launch.type,'LH'),
0092             iw_list = [iw_list, iw];
0093             %
0094             if ~isfield(waves{iw},'a_sdNpar'),
0095                 disp(['a_sdNpar not specified is set to zero for wave ',int2str(iw)])
0096                 a_sdNpar(iw) = 0;% used only in abslim loop
0097             else
0098                 a_sdNpar(iw) = waves{iw}.a_sdNpar;
0099             end
0100         end
0101     end
0102     %
0103     for ik = 1:nk_abs,
0104         %
0105         for iw = iw_list;            
0106             %
0107             for ib = 1:nb(iw),
0108                 waves{iw}.rays{ib}.sdNpar = waves{iw}.rays{ib}.sdNpar(1) +ik*a_sdNpar(iw)*waves{iw}.rays{ib}.ss;
0109             end
0110             %
0111             if isfield(waves{iw},'waveparam') && isfield(waves{iw}.waveparam,'n_rf_list'),
0112                 n_rf_list = waves{iw}.waveparam.n_rf_list;
0113             elseif isfield(dkeparam,'n_rf_list'),
0114                 n_rf_list = dkeparam.n_rf_list;
0115             else
0116                 disp('Warning : no n_rf_list data found, enforced to 0:3')
0117                 n_rf_list = 0:3;
0118             end
0119             %
0120             % the wave linear deposition is recalculated accounting for a_sdNpar effect
0121             %
0122             waves{iw}.a_sdNpar = ik*a_sdNpar(iw);
0123             %
0124             waves{iw} = wave_process_jd(equil,waves{iw},n_rf_list);%Wave processing and relativistic linear damping
0125             %
0126         end
0127         %
0128         % the radial grid is recalculated accounting for a_sdNpar effect
0129         %
0130         [dkeparam.rho_S,dkeparam.pnmax_S] = optgrid_dke_jd(waves,rho_S_init,dkeparam.pnmax_S,equil,w_mask,dkedisplay.display_mode);
0131         %
0132         [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0133         %
0134         ny = size(ZP0(end).xy_rf_fsav,2);
0135         %
0136         yp_rf = 2*pi*equil.Rp*sum(ZP0(end).xy_rf_fsav*mksa.P_ref.*repmat(equilDKE.xdV_2piRp_dke.',[1,ny]),1);
0137         %
0138         P0 = 0;
0139         Pabs = 0;
0140         yb = dke_out.yb;
0141         %
0142         for iw = iw_list,
0143             %
0144             ibmin = sum(nb(1:iw-1)) + 1;
0145             ibmax = sum(nb(1:iw)) + 1;
0146             %
0147             if iw < nw
0148                 ymask = yb(ibmin):(yb(ibmax) - 1);
0149             else
0150                 ymask = yb(ibmin):ny;
0151             end            
0152             %
0153             Pabs = Pabs + 1e6*sum(yp_rf(ymask));
0154             for ib = 1:nb(iw),
0155                 P0 = P0 + waves{iw}.rays{ib}.P0_2piRp*2*pi*equil.Rp;
0156             end
0157         end
0158         %
0159         if P0 > 0 && isfield(dkeparam,'abs_lim') && Pabs/P0 < dkeparam.abs_lim,
0160             %
0161             if a_sdNpar(iw) == 0,
0162                 %
0163                 a_sdNpar(iw) = 0.01;
0164                 %
0165             end
0166             %
0167         else
0168             break
0169         end
0170     end
0171 % *************************************************************************
0172 %
0173 % Case with sdNpar increasing at each caustic point
0174 %
0175 elseif dkeparam.rt_mode == 5,
0176     %
0177     for iw = 1:nw,
0178         for ib = 1:nb(iw),
0179             scaustic = logical([0,diff(waves{iw}.rays{ib}.pmode) ~= 0]);
0180             sscattering = cumsum(scaustic);
0181             %
0182             waves{iw}.rays{ib}.sdNpar = waves{iw}.rays{ib}.sdNpar(1)*waves{iw}.b_sdNpar.^sscattering;
0183         end
0184     end
0185     %
0186     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0187 %
0188 % *************************************************************************
0189 %
0190 % Seek mode : optimize according to highest driven current
0191 %
0192 elseif dkeparam.rt_mode == 3,%
0193     %
0194     if isfield(dkeparam,'seekparam'),
0195         %
0196         submode = dkeparam.seekparam.submode;
0197         %
0198         if submode == 1,%case where sdNpar/sNar is a constant factor
0199             val = dkeparam.seekparam.f_sdNpar_in;%This is in dkeparam because it is the same for all waves
0200             step = dkeparam.seekparam.f_step;
0201         elseif submode == 2,
0202             val = dkeparam.seekparam.a_sdNpar_in;%This is in dkeparam because it is the same for all waves
0203             step = dkeparam.seekparam.a_step;
0204             %
0205             sdNpar_in = cell(1,nw);
0206             for iw = 1:nw,
0207                 sdNpar_in{iw} = zeros(1,nb(iw));
0208                 for ib = 1:nb(iw),
0209                     sdNpar_in{iw}(ib) = waves{iw}.rays{ib}.sdNpar(1);
0210                 end
0211             end
0212         else
0213             error('unknown method in rt_mode 3')
0214         end
0215     else
0216         error('dkeparam.seekparam must be specified for rt_mode 3')
0217     end
0218     %
0219     opt_loop = 0;
0220     eta = 0;
0221     firststep = 1;
0222     %
0223     while opt_loop == 0,
0224         %
0225         eta_old = eta;
0226         %
0227         for iw = 1:nw,
0228             for ib = 1:nb(iw),
0229                 if submode == 1,
0230                     waves{iw}.rays{ib}.sdNpar = abs(real(waves{iw}.rays{ib}.sNpar))*val;
0231                 elseif submode == 2,
0232                     waves{iw}.rays{ib}.sdNpar = sdNpar_in{iw}(ib) + waves{iw}.rays{ib}.ss*val;
0233                 end
0234             end
0235         end
0236         %
0237         dNstr = num2str(val);
0238         dNstr = dNstr(1:min(5,length(dNstr)));
0239         %
0240         for is=1:length(dNstr),if dNstr(is) == '.',dNstr(is) = 'p';end,end
0241         %
0242         % eval(['diary ',path_simul,'DKE_DIARY_',equil.id,'_',simul.id,'_',dNstr,'.txt']);
0243         %
0244         [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0245         %
0246         % CD efficiency (A/W)
0247         %
0248         Prf = sum(ZP0.x_rf_fsav(:).*equilDKE.xdV_2piRp_dke(:)*2*pi*equilDKE.Rp*mksa.P_ref);
0249         Itot = sum(Zcurr.x_tot_fsav(:).*equilDKE.xdA_dke(:)*mksa.j_ref);
0250         eta = abs(Itot/Prf);
0251         %
0252         eval(['save ',path_simul,'DKE_RESULTS_',equil.id,'_',simul.id,'_',dNstr,'.mat']);
0253         %
0254         disp(['--> Step ',dNstr,' calculated. eta = ',num2str(eta),' A/W'])
0255         %
0256         % diary off
0257         %
0258         % This function calculates the power deposited along each ray using LUKE results.
0259         %
0260         waves = alphaphi_luke_jd(waves,dke_out,gridDKE);
0261         %
0262         if eta > eta_old,
0263             %
0264             val = val + step;
0265             %
0266             if eta_old > 0,
0267                 firststep = 0;
0268             end
0269             %
0270         elseif firststep == 1,
0271             %
0272             step = - step;
0273             val = val + 2*step;
0274             %
0275             firststep = 0;
0276             %
0277         else
0278             %
0279             opt_loop = 1;
0280             %
0281             val = val - step;
0282             %
0283         end
0284         %
0285         if opt_loop == 0 && val < -1000*eps,
0286             %
0287             opt_loop = 1;
0288             %
0289             val = val - step;
0290             %
0291         end
0292     end
0293     %
0294     % eval(['diary ',path_simul,'DKE_DIARY_',equil.id,'_',simul.id,'.txt']);
0295     %
0296     for iw = 1:nw,
0297         for ib = 1:nb(iw),
0298             if submode == 1,
0299                 waves{iw}.rays{ib}.sdNpar = abs(waves{iw}.rays{ib}.sNpar)*val;
0300             elseif submode == 2,
0301                 waves{iw}.rays{ib}.sdNpar = sdNpar_in{iw}(ib) + waves{iw}.rays{ib}.ss*val;
0302             end
0303         end
0304     end
0305     %
0306     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct);
0307     %
0308     % diary off
0309 elseif dkeparam.rt_mode == -1,
0310     %
0311     % No FP calculations, return equil and waves
0312     %
0313     Znorm = [];
0314     Zcurr.x_0_fsav = [];
0315     ZP0.x_rf_fsav = [];
0316     %
0317     dke_out.equil = equil;
0318     dke_out.waves = waves;
0319     dke_out.simul = simul;
0320     %
0321     dke_out.XXf0 = {[]};
0322     dke_out.XXfM = [];
0323     dke_out.ZXXD_a = [];
0324     dke_out.ZXXF_a = [];
0325     dke_out.ZXXD_c = [];
0326     dke_out.ZXXF_c = [];
0327     dke_out.ZXXD_e = [];
0328     dke_out.ZXXF_e = [];
0329     dke_out.XXILor = [];
0330     dke_out.xepsi = [];
0331     dke_out.xEfield_validity = [];
0332     dke_out.FPtime = [];
0333     %
0334     radialDKE = [];
0335     equilDKE = [];
0336     momentumDKE = [];
0337     gridDKE = [];
0338     Zmomcoef = [];
0339     Zbouncecoef = [];
0340     Zmripple = [];
0341     mksa = [];
0342     XXsinksource = [];  
0343     dke_warnings = [];
0344     %
0345 else
0346     Znorm = [];
0347     Zcurr = [];
0348     ZP0 = [];
0349     dke_out = [];
0350     radialDKE = [];
0351     equilDKE = [];
0352     momentumDKE = [];
0353     gridDKE = [];
0354     Zmomcoef = [];
0355     Zbouncecoef = [];
0356     Zmripple = [];
0357     mksa = [];
0358     XXsinksource = [];  
0359     dke_warnings = [];
0360     warning('rt mode not recognized in loop_run_lukert.m');
0361 end
0362 %
0363 
0364 
0365 
0366 
0367 
0368 
0369 
0370 
0371 
0372 
0373 
0374 
0375 
0376 
0377 
0378 
0379 
0380