Description of cronos2luke

0001 function [sortie,memoire] = cronos2luke(par,equi,neo,gene,prof,compo,impur,proto,memoire)
0002 %
0003 %  Non-thermal distribution and associated current (LH, EC, runaways, ...) calculated  by the codes C3PO/LUKE for CRONOS
0004 %
0005 % INPUTS :
0006 %
0007 %   - par     : parametres C3PO/LUKE
0008 %   - equi    : donnees de l'equilibre plasma (data.equi)
0009 %   - neo     : donnees neoclassiques (data.neo)
0010 %   - gene    : parametres generaux (param.gene)
0011 %   - prof    : profils des donnees calculees par le code (data.prof)
0012 %   - compo   : composition du plasma: charge et masse des atomes (param.compo)
0013 %   - impur   : sous strcuture des impurtes (data.impur)
0014 %   - proto   : prototype de la structure pour les sources, valeurs a zeros (proto : zsourceproto;)
0015 %   - memoire : structure des dernieres valeurs calculees
0016 %
0017 % OUTPUTS :
0018 %
0019 %     sortie  : structure de type source remplie par le module (sortie = proto)
0020 %     memoire : datak.memoire.hyb (valeur de reference pour le dernier calcul complet, ne pas utiliser dans cette fonction, reserve pour d'autres modules)
0021 %
0022 % by J. Decker, Y. Peysson, J.-F. Artaud and V. Basiuk
0023 %
0024 sortie     = proto;% output initialization
0025 sortie.err = 0;% flag d'erreur a 0 par defaut (pas d'erreur)
0026 %
0027 %--------------------------- LIEN AVEC C3PO + LUKE ------------------
0028 %
0029 % tokamak name
0030 %
0031 if ~isfield(par,'tok')
0032     tokamak = 'unknown';
0033 else
0034     tokamak = par.tok;
0035 end
0036 %
0037 % shot number
0038 %
0039 if ~isfield(par,'shot')
0040     shotnum = '0';
0041 else
0042     shotnum = num2str(par.shot); 
0043 end
0044 %
0045 % shot time
0046 %
0047 if isempty(memoire.data),
0048     shotime = par.temps_cons.tempsini;% the distribution is evolved between -infinity and the initial time
0049 else
0050     shotime = par.temps_cons.tempsend;% the distribution is evolved between the initial and the final time
0051 end
0052 %
0053 shotimestr = num2str(shotime,'%6.4f');
0054 %
0055 % equil ID
0056 %
0057 equil_id = [tokamak,'_',shotnum,'_',shotimestr];
0058 %
0059 % PARAMETRES LUKE
0060 %
0061 dkeparam = make_dkeparam_cronos(par);
0062 %
0063 % PATH LUKE
0064 %
0065 dkepath = load_structures_yp('dkepath','','');
0066 %
0067 % PARAMETRES DISPLAY
0068 %
0069 dkedisplay = make_display_cronos;
0070 dkedisplay.rho_display  = par.rho_display;%Approximate radial position at which results are displayed if display.display_mode >= 1 [1,1]. Not used if NaN
0071 dkedisplay.display_mode = par.display_mode;%Display mode: (0): no, (1): partial, (2): full
0072 %
0073 % PARAMETRES EQUILIBRE
0074 %
0075 equil = make_equil_cronos(equil_id,equi,gene,prof,compo,impur);
0076 %
0077 % PARAMETRES OHMIC, transformation en flux poloidal puis en coordonnees radiales de LUKE
0078 %
0079 % psipncronos    : coordonnee radiale normalisee de CRONOS en flux poloidal
0080 % psipncronoshel : coordonnee radiale normalisee en flux poloidal associee a HELENA pour BR BZ et BPHI
0081 % psitncronoshel : coordonnee radiale normalisee en flux toroidal associee a HELENA pour BR BZ et BPHI
0082 % psipn_jd       : coordonnee radiale normalisee en flux poloidal associee a LUKE
0083 %
0084 nrho_ohm       = 101;
0085 ohm.rho        = linspace(0,1,nrho_ohm);% grille radiale rhoG de LUKE
0086 %psipncronos    = abs(equi.psi-equi.psi(1))/max(abs(equi.psi-equi.psi(1)));
0087 %
0088 if equi.psiRZ(1) == equi.psiRZ(2)
0089    disp('non monotonic psiRZ');
0090    psiRZ = interp1(gene.x*equi.rhomax,equi.psi,double(equi.rhoRZ),'pchip','extrap');
0091    equi.psiRZ = double(psiRZ);
0092 else
0093    psiRZ = double(equi.psiRZ); 
0094 end
0095 %
0096 psipncronoshel = abs(psiRZ-psiRZ(1))/max(abs(psiRZ-psiRZ(1)));
0097 psitncronoshel = equi.rhoRZ/equi.rhomax;
0098 psipn_ohm       = rho2psi_jd(equil,ohm.rho);
0099 %---------------------------------------------------------
0100 % calcul du champ // pour LUKE (avec gestion des signes)
0101 % passage de prof.epar de la grille en x sur la grille HELENA
0102 % puis passage sur la grille LUKE en faisant correspondre rhohel et psihel
0103 % puis multiplication par le facteur geometrique (pour s'assurer que l'integrale redonne le courant)
0104 %
0105 % prof.epar (V/m): signe par rapport à IP
0106 % ohm.epsi (V/m): signe par rapport à BT
0107 % - signe is electron charge
0108 %----------------------------------------------------------
0109 epsicronoshel     = -sign(gene.signe.ip)*sign(gene.signe.b0)*interp1(gene.x,prof.epar,psitncronoshel);% electric field with respect to BT on helena grid
0110 ohm.epsi     = complex(interp1(psipncronoshel,epsicronoshel,psipn_ohm));%Imaginary for considering the definition of the flux surface averaging of CRONOS in LUKE
0111 %
0112 % PARAMETRES RIPPLE
0113 %
0114 ripple.mode = par.ripple;
0115 %
0116 % PARAMETRES FAST ELECTRON TRANSPORT
0117 %
0118 if par.Dr0 == 0 && par.Vr0 == 0
0119     %
0120     % cas sans transport, structure vide
0121     %
0122     transpfaste = '';
0123 else
0124     transpfaste.vparmin  = par.vparmin_r;%Lower limit of the parallel velocity dependence of the radial diffusion and pinch (vth_ref or vth)
0125     transpfaste.norm_ref = par.norm_ref;%Normalization procedure for the lower limit of radial transport parallel velocity dependence: (0) from local values Te and ne, (1) from reference values Te_ref and ne_ref
0126     transpfaste.Dr0      = par.Dr0;%Core radial diffusion coefficient (m^2/s)
0127     transpfaste.pDr      = par.pDr;%Coefficient for the Dr radial profile
0128     transpfaste.Vr0      = par.Vr0;%Core radial pinch coefficient (m/s)
0129     transpfaste.pVr      = par.pVr;%Coefficient for the Vr radial profile
0130     transpfaste.Dr_model = par.Dr_model;%Radial diffusion type: (0) no vpar dependence, (1) magnetic turbulence model (vpar/vth),
0131     transpfaste.Vr_model = par.Vr_model;%Radial pinch type: (0) no vpar dependence, (1) magnetic turbulence model (vpar/vth),
0132 end
0133 %
0134 % PARAMETRES WAVE LH + EC (C3PO)
0135 %
0136 if ~isfield(par,'conf')
0137  par.conf = tokamak;
0138 end 
0139 %
0140 cons=par.cons;
0141 if strcmp(par.luke_mode,'lh')
0142   [launchs_lh] = make_wavelaunchs_cronos_lh(gene,par,cons.lh_cons);%for the LH waves
0143   launchs_ec   = [];
0144 end
0145 if strcmp(par.luke_mode,'ec')
0146   [launchs_ec] = make_wavelaunchs_cronos_ec(equil,par,cons.ecrh_cons);%for the EC waves
0147   launchs_lh   = [];
0148 end
0149 if strcmp(par.luke_mode,'mixte')
0150   [launchs_lh] = make_wavelaunchs_cronos_lh(gene,par,cons.lh_cons);%for the LH waves
0151   [launchs_ec] = make_wavelaunchs_cronos_ec(equil,par,cons.ecrh_cons);%for the EC waves
0152 end
0153 %
0154 launchs = [launchs_lh,launchs_ec];
0155 %
0156 mask = [];
0157 prf = 0;
0158 for iw = 1:length(launchs),
0159     %
0160     if strcmp(launchs{iw}.type,'LH')
0161         wprf = sum(launchs{iw}.bPlhtot);
0162     elseif strcmp(launchs{iw}.type,'EC')
0163         wprf = sum(launchs{iw}.yP_L);
0164     end
0165     %
0166     if ~isfield(par,'prfmin') || wprf >= par.prfmin,
0167         mask = [mask,iw];
0168         prf = prf + wprf;
0169     else
0170         disp(['wave #',int2str(iw),': dismissed as Prf = ',num2str(wprf),' < Prfmin = ',num2str(par.prfmin),'.'])
0171     end
0172 end
0173 launchs = launchs(mask);  
0174 nw = length(launchs);
0175 %
0176 if nw == 0
0177    sortie.j                    = zeros(size(gene.x));
0178    sortie.el                   = zeros(size(gene.x));
0179    memoire.data.output.j       = sortie.j;
0180    memoire.data.output.el      = sortie.el;
0181    memoire.data.output.paddrf  = 0;
0182    memoire.data.output.prf     = 0;
0183    memoire.data.output.ne      = prof.ne;
0184    memoire.data.output.irf     = 0;
0185    return 
0186 end
0187 wavestructs = cell(1,length(launchs));
0188 %
0189 dkeparam.w_mask = [];% waves to be accounted for in radial grid optimization: EC only or LH if no EC present
0190 %
0191 for iw = 1:nw,
0192     wavestructs{iw} = make_wavert_cronos(launchs{iw},par);%parameter for ray-tracing are adapted to the selected wave
0193     %
0194     wavestructs{iw}.equil_id = equil.id;
0195     %
0196     if strcmp(launchs{iw}.type,'EC')
0197         %
0198         dkeparam.w_mask = [dkeparam.w_mask,iw];
0199         disp(['wave #',int2str(iw),': type :',launchs{iw}.type])
0200         %
0201     else
0202         disp(['wave #',int2str(iw),': type :',launchs{iw}.type,' n// = ',num2str(wavestructs{iw}.launch.bNpar0,4)])
0203         disp(['wave #',int2str(iw),': type :',launchs{iw}.type,' Plh = ',num2str(wavestructs{iw}.launch.bPlhtot,4)])
0204     end
0205     %
0206 end
0207 if isempty(dkeparam.w_mask)
0208     dkeparam.w_mask = 1:length(launchs);
0209 end
0210 %
0211 % accessibility test
0212 %
0213 if nw > 0
0214   acces_lh
0215 end
0216 %
0217 distrib_RF  = [];
0218 distrib_TOT = [];
0219 distrib_OHM = [];
0220 %
0221 if ~isempty(memoire)
0222   if isfield(memoire,'data')
0223     if isfield(memoire.data,'edistribution')  
0224         disp(['-------------using the last distribution function------------'])
0225         tm = memoire.data.edistribution(end).shotime;
0226         if tm == shotime,% remove last distribution info if time has not advanced
0227             memoire.data.edistribution = memoire.data.edistribution(1:end-1);
0228         end
0229         %
0230         if length(memoire.data.edistribution) >= 1,% the time here is necessarily anterior
0231             distrib_RF     = memoire.data.edistribution(end).f0_RF;
0232             distrib_TOT    = memoire.data.edistribution(end).f0_TOT;
0233             distrib_OHM    = memoire.data.edistribution(end).f0_OHM;
0234             dkeparam.rho_S = memoire.data.edistribution(end).rho_S;
0235             dkeparam.tn    = (shotime - memoire.data.edistribution(end).shotime)*distrib_RF.mksa.nhu_ref;% LUKE time step
0236         end
0237     end
0238   end
0239 end
0240 %
0241 if par.save == 2,%save before LUKE
0242    %
0243   try
0244     workdir = fileparts(evalin('base','param.gene.file'));
0245   catch 
0246     workdir = [];
0247   end
0248   if isempty(workdir)
0249     disp('you can not save the luke file in batch mode')
0250   else
0251     save([workdir,'/before_luke@',equil_id,'.mat']);
0252     error('Exiting after saving LUKE parameters')
0253   end
0254 
0255 end
0256 %
0257 %--------------------------- Execution du trace de rayon C3PO et du code cinetique LUKE ---------------------------
0258 %
0259 ik = 1;
0260 warning off
0261 disp(['LUKE is running ....'])
0262 %
0263 signe_out = -gene.signe.ip;
0264 %
0265 % CALCULATION WITH RF ONLY
0266 %
0267 [Znorm_RF,Zcurr_RF,ZP0_RF,dke_out_RF,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = ...
0268     run_lukert_batchopt('cronos',wavestructs,[],dkepath,equil,dkeparam,dkedisplay,[],transpfaste,ripple,distrib_RF,par);
0269 %
0270 j_RF_c = signe_out*Zcurr_RF.x_0_vcfsav*mksa.j_ref*1e6;
0271 I_RF_c = sum(signe_out*Zcurr_RF.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref*1e6);
0272 p_RF_c = ZP0_RF.x_rf_fsav.'*mksa.P_ref*1e6;
0273 P_RF_c = sum(p_RF_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0274 f0_RF = dke_out_RF.Zf0_interp;%For time evolution
0275 %
0276 waves = dke_out_RF.waves;%waves reused for synergy calculations
0277 dkeparam.rho_S = radialDKE.xrho_S_dke;% retrieve RF rho grid if (re)calculated within run_lukert
0278 dkeparam.rt_mode = 1;%a_sdNpar fixed by previous RF calculation
0279 %
0280 if strcmp(par.syn,'Yes') && sum(abs(ohm.epsi)) ~= 0,%Calculate the synergy between RF and OHMIC
0281     %
0282     % CALCULATION WITH RF AND OHMIC
0283     %
0284     [Znorm_TOT,Zcurr_TOT,ZP0_TOT,dke_out_TOT] = run_lukert_batchopt('cronos',[],waves,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,distrib_TOT,par);
0285     %
0286     j_TOT_c = signe_out*Zcurr_TOT.x_0_vcfsav*mksa.j_ref*1e6;
0287     I_TOT_c = sum(signe_out*Zcurr_TOT.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref*1e6);
0288     p_RF_TOT_c = ZP0_TOT.x_rf_fsav.'*mksa.P_ref*1e6;
0289     % p_OHM_TOT_c = ZP0_TOT.x_e_fsav*mksa.P_ref*1e6;
0290     % P_RF_TOT_c = sum(p_RF_TOT_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0291     % P_OHM_TOT_c = sum(p_OHM_TOT_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0292     %
0293     f0_TOT = dke_out_TOT.Zf0_interp;%For time evolution
0294     %
0295     % CALCULATION OHMIC ONLY
0296     %
0297     [Znorm_OHM,Zcurr_OHM,ZP0_OHM,dke_out_OHM] = run_lukert_batchopt('cronos',[],[],dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,distrib_OHM,par);
0298     %
0299     j_OHM_c = signe_out*Zcurr_OHM.x_0_vcfsav*mksa.j_ref*1e6;
0300     I_OHM_c = sum(signe_out*Zcurr_OHM.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref*1e6);
0301     % p_OHM_c = ZP0_OHM.x_e_fsav*mksa.P_ref*1e6;
0302     % P_OHM_c = sum(p_OHM_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0303     %
0304     f0_OHM = dke_out_OHM.Zf0_interp;%For time evolution
0305     %
0306 else
0307     %
0308     j_TOT_c     = j_RF_c;% total current is RF current
0309     I_TOT_c     = I_RF_c;
0310     p_RF_TOT_c  = p_RF_c;% RF deposition profile in the absence of ohmic
0311     j_OHM_c     = zeros(1,length(j_RF_c));% no ohmic current
0312     I_OHM_c     = 0;
0313     f0_TOT      = f0_RF;% electron distribution in the absence of ohmic
0314     f0_OHM      = [];% maxwellian distribution
0315     %
0316 end
0317 %
0318 % SYNERGISTIC EFFECTS - proportional to the electric field in linear limit
0319 %
0320 % Note : p_OHM_SYN_c is commented as it should be already accounted for in cronos E.J with J including j_SYN_c
0321 %
0322 j_SYN_c = j_TOT_c - j_RF_c - j_OHM_c;% synergistic current
0323 I_SYN_c = I_TOT_c - I_RF_c - I_OHM_c;% synergistic current
0324 % p_OHM_SYN_c = p_OHM_TOT_c - p_OHM_c;% synergistic ohmic heating
0325 %
0326 warning on
0327 %
0328 if strcmp(dkedisplay.display_mode,'PARTIAL_VISUAL') || strcmp(dkedisplay.display_mode,'VISUAL'),
0329     %
0330     visu_luke
0331     %
0332 end
0333 %
0334 % par.ipfrac : maximum fraction of IP allowed for RF current
0335 %
0336 if isfield(par,'ipfrac'),
0337     %
0338     I_boot = zintsurf(neo.jboot,gene.x,equi.spr,equi.rhomax);
0339     %
0340     if (I_RF_c+I_boot)/equi.ip > par.ipfrac,
0341         %
0342         normip = par.ipfrac*equi.ip/(I_RF_c+I_boot);
0343         disp(['----> to avoid I_L_H + I_b_o_o_t > Ip'])
0344         disp(['normalisation =',num2str(normip,3)])
0345         %
0346     else
0347         %
0348         normip = 1;
0349         %
0350     end
0351 end
0352 %
0353 sortie.err = normip;
0354 %
0355 %--------------------------- REMPLISSAGE DE LA STRUCTURE DE SORTIE POUR CRONOS ---------------------------
0356 %
0357 %
0358 % preparation des variables de sortie (sortie.j, sortie.el) pour CRONOS
0359 %
0360 method = 'linear';
0361 %
0362 xrhoT = [0,equilDKE.xrhoT,1];
0363 j_RF_c = [j_RF_c(1),j_RF_c,0];
0364 j_SYN_c = [j_SYN_c(1),j_SYN_c,0];
0365 p_RF_TOT_c = [p_RF_TOT_c(1),p_RF_c,0];
0366 %
0367 cj_RF_c  = interp1(xrhoT,j_RF_c,gene.x,method);% (A/m2)
0368 cj_SYN_c  = interp1(xrhoT,j_SYN_c,gene.x,method);% (A/m2)
0369 cp_RF_TOT_c = interp1(xrhoT,p_RF_TOT_c,gene.x,method);% (W/m3)
0370 %
0371 % remise de la puissance et du courant de LUKE
0372 %
0373 irfnew = zintsurf(cj_RF_c,gene.x,equi.spr,equi.rhomax);
0374 isynnew = zintsurf(cj_SYN_c,gene.x,equi.spr,equi.rhomax);
0375 paddrfnew = zintvol(cp_RF_TOT_c,gene.x,equi.vpr,equi.rhomax);
0376 %
0377 if abs(irfnew - I_RF_c)/abs(I_RF_c) > 0.01,
0378     disp(['WARNING : surface integral on LUKE and CRONOS grids give different results. RF current renormalized by ',num2str(I_RF_c / irfnew)])
0379 end
0380 if abs(isynnew - I_SYN_c)/max([abs(I_SYN_c),eps]) > 0.01,
0381     disp(['WARNING : surface integral on LUKE and CRONOS grids give different results. SYN current renormalized by ',num2str(I_SYN_c / isynnew)])
0382 end
0383 if abs(paddrfnew - P_RF_c)/abs(P_RF_c) > 0.01,
0384     disp(['WARNING : volume integral on LUKE and CRONOS grids give different results. RF power renormalized by ',num2str(P_RF_c / paddrfnew)])
0385 end
0386 cj_RF_c  = cj_RF_c  * I_RF_c / irfnew * normip;
0387 if isynnew ~= 0,
0388     cj_SYN_c  = cj_SYN_c  * I_SYN_c / isynnew;
0389 end
0390 if paddrfnew ~= 0,
0391     cp_RF_TOT_c = cp_RF_TOT_c * P_RF_c / paddrfnew;
0392 end
0393 %
0394 % SMOOTHING
0395 %
0396 if isfield(par,'smooth') && strcmp(par.smooth,'Yes')
0397     newj  = smooth_dke_yp(cj_RF_c.',par.width_smooth,par.method_smooth).';  % diffusion phenomenologique du courant
0398     norm1 = zintsurf(cj_RF_c,gene.x,equi.spr,equi.rhomax);
0399     norm2 = zintsurf(newj,gene.x,equi.spr,equi.rhomax);
0400     cj_RF_c = newj/norm2*norm1;
0401     %
0402     newjsyn  = smooth_dke_yp(cj_SYN_c.',par.width_smooth,par.method_smooth).';  % diffusion phenomenologique du courant
0403     norm1 = zintsurf(cj_SYN_c,gene.x,equi.spr,equi.rhomax);
0404     norm2 = zintsurf(newjsyn,gene.x,equi.spr,equi.rhomax);
0405     if norm1 ~= 0,
0406         cj_SYN_c  = newjsyn/norm2*norm1;
0407     end
0408     %
0409     newel  = smooth_dke_yp(cp_RF_TOT_c.',par.width_smooth,par.method_smooth).';  % diffusion phenomenologique du courant
0410     norm1 = zintvol(cp_RF_TOT_c,gene.x,equi.vpr,equi.rhomax);
0411     norm2 = zintvol(newel,gene.x,equi.vpr,equi.rhomax);
0412     if norm1 ~= 0,
0413         cp_RF_TOT_c = newel/norm2*norm1;
0414     end
0415 end
0416 %
0417 % if negative RF power is found at any xy grid point, it means the calculation was somewhat unstable
0418 %
0419 if any(any(ZP0_RF.xy_rf_fsav < 0)),
0420     %
0421     sortie.err = sortie.err + 1i;
0422     %
0423     disp('Problem in LUKE, negative absorbed power calculated for some value, use the previous LUKE calculation if possible')
0424     %
0425     if isfield(memoire.data,'output')
0426         %
0427         % assume constant efficiency n.I/P
0428         %
0429         cj_RF_c   = memoire.data.output.j*prf/memoire.data.output.prf.*memoire.data.output.ne./prof.ne;
0430         cj_SYN_c = 0;
0431         cp_RF_TOT_c  = memoire.data.output.el*prf/memoire.data.output.prf;
0432         %
0433     end
0434 end    
0435 %
0436 % In principle the synergy increases with both RF power and efield. Thus it
0437 % makes sense, to avoid divergence in sigmasyn, to cap the sigmasyn with
0438 % its value obtained at the max of abs(cj_SYN_c)
0439 %
0440 sigmasyn = cj_SYN_c./prof.epar;
0441 sigmasyn(prof.epar == 0) = 0;% To avoid NAN
0442 sigmasynmax = sigmasyn(find(abs(cj_SYN_c) == max(abs(cj_SYN_c)),1,'first'));%value obtained at the max of abs(cj_SYN_c)
0443 sigmasyn = sign(sigmasyn).*min([abs(sigmasyn);repmat(abs(sigmasynmax),[1,length(sigmasyn)])]);% caping
0444 %
0445 sortie.j = cj_RF_c;% RF current in the absence of ohmic
0446 sortie.sigmasyn = sigmasyn;% synergistic conductivity (= HOT - COLD).
0447 sortie.el = cp_RF_TOT_c;% Note : the RF power deposition profile calculated in the presence of ohmic must be used.
0448 %
0449 irf = zintsurf(cj_RF_c,gene.x,equi.spr,equi.rhomax);
0450 paddrf = zintvol(cp_RF_TOT_c,gene.x,equi.vpr,equi.rhomax);
0451 %
0452 disp(['I_L_H = ',num2str(irf/1e3,3),' kA / I_p =',num2str(equi.ip/1e3,4),' kA'])
0453 disp(['P_L_H = ',num2str(paddrf/1e3,3),' kW / P_r_e_f = ',num2str(prf/1e3,3),' kW'])
0454 %
0455 memoire.data.output.j         = cj_RF_c;
0456 memoire.data.output.el        = cp_RF_TOT_c;
0457 memoire.data.output.paddrf    = paddrf;
0458 memoire.data.output.prf       = prf;
0459 memoire.data.output.ne        = prof.ne;
0460 memoire.data.output.irf       = irf;
0461 %
0462 for iw = 1:nw,
0463     if isfield(wavestructs{iw}.launch,'aloha') & wavestructs{iw}.launch.aloha.mode ~= 0
0464        if isfield(memoire.data.output,'npar') 
0465            if iw <= length(memoire.data.output.npar)
0466              memoire.data.output.npar{iw}(end+1,:) = wavestructs{iw}.launch.bNpar0;
0467            else
0468              memoire.data.output.npar{iw}(1,:) = wavestructs{iw}.launch.bNpar0;
0469            end
0470        else
0471            memoire.data.output.npar{iw}(1,:) = wavestructs{iw}.launch.bNpar0;
0472        end
0473        if isfield(memoire.data.output,'pow') 
0474            if iw <= length(memoire.data.output.pow)
0475              memoire.data.output.pow{iw}(end+1,:) = wavestructs{iw}.launch.bPlhtot;
0476            else
0477              memoire.data.output.pow{iw}(1,:) = wavestructs{iw}.launch.bPlhtot;
0478            end
0479        else
0480            memoire.data.output.pow{iw}(1,:) = wavestructs{iw}.launch.bPlhtot;
0481        end
0482     end
0483 end
0484 %
0485 if isfield(memoire.data,'edistribution')
0486     memoire.data.edistribution(end + 1).shotime = shotime;%For time evolution
0487 else
0488     memoire.data.edistribution.shotime = shotime;%For time evolution
0489 end
0490 %
0491 memoire.data.edistribution(end).f0_RF   = f0_RF;%For time evolution
0492 memoire.data.edistribution(end).f0_TOT  = f0_TOT;%For time evolution
0493 memoire.data.edistribution(end).f0_OHM  = f0_OHM;%For time evolution
0494 memoire.data.edistribution(end).rho_S   = radialDKE.xrho_S_dke;
0495 %
0496 %
0497 % For diagnostic post-processing (HXR,...)
0498 %
0499 %memoire.data.edistribution.dkeparam  = dkeparam;%For HXR calculations
0500 %memoire.data.edistribution.dkedisplay = momentumDKE.dkedisplay;%For HXR calculations
0501 %memoire.data.edistribution.equil = equil;%For HXR calculations
0502 %memoire.data.edistribution.radialDKE = radialDKE;%For HXR calculations
0503 %memoire.data.edistribution.Zcurr = Zcurr;%For HXR calculations
0504 %memoire.data.edistribution.hxr = hxr;%For HXR calculations
0505 %memoire.data.edistribution.hxrparam = hxrparam;%For HXR calculations
0506 %memoire.data.edistribution.mksa = mksa;%For HXR calculations
0507 %memoire.data.edistribution.momentumDKE = momentumDKE;%For HXR calculations
0508 %memoire.data.edistribution.dke_out = dke_out_RF;%For HXR calculations
0509 %
0510 if par.save == 1,
0511   try
0512     workdir = fileparts(evalin('base','param.gene.file'));
0513   catch 
0514     workdir = [];
0515   end
0516   if isempty(workdir)
0517     disp('you can not save the luke file in batch mode')
0518   else
0519     save([workdir,'/last_luke@',equil_id,'.mat']);
0520   end
0521 end
0522 end
0523 %
0524 % Generic function to use either direct call to run_lukert or batch run
0525 %
0526 function [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = ...
0527     run_lukert_batchopt(id_simul,wavestructs,waves,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,distrib,par)
0528     %
0529     dkepath.remnum = par.remnum;% remote profile number. Use 0 for local sequential calculations
0530     dkepath.remtimout = par.remtimout;% timeout for remote calculations, in minutes. Use (0) to return later and check results by running this script again
0531     dkepath.rempause = par.rempause;% pause between job checks (in minutes) (only for remtimout > 0)
0532     dkepath.debugmode = par.debugmode;% debug mode (0) off (1) on
0533     dkepath.clean = par.clean;% clean job files/folders : both on remote and local (2), only on remote (1) or not at all (0)
0534     %
0535     % backup & save options
0536     %
0537     opt.save = par.saveluke;% (1): save LUKE_RESULTS or return error (0): return results or error
0538     if isinf(par.fields),
0539         opt.fields = NaN;% selected fields in returned lukestructs. Use NaN for all fields
0540     else
0541         opt.fields = par.fields;% selected fields in returned lukestructs. Use NaN for all fields
0542     end
0543     opt.backup = par.backup;% (1) save fluctuation time steps (0) do not save
0544     opt.waves = par.waves;% (1) save waves at each time step (0) do not save
0545     %
0546     lukestructs = struct;
0547     lukestructs.simul.id_simul = id_simul;
0548     lukestructs.simul.path = '';
0549     lukestructs.dkeparam = dkeparam;
0550     lukestructs.dkepath = dkepath;
0551     lukestructs.dkedisplay = dkedisplay;
0552     lukestructs.equil = equil;
0553     lukestructs.waves = waves;
0554     lukestructs.wavestructs = wavestructs;
0555     lukestructs.ohm = ohm;
0556     lukestructs.transpfaste = transpfaste;
0557     lukestructs.ripple = ripple;
0558     lukestructs.Zf0_interp = distrib;
0559     lukestructs.opt = opt;
0560     %
0561     lukestructs = run_lukert(lukestructs,dkepath);
0562     %
0563     if isfield(lukestructs,'err'),
0564         zverbose('Remote LUKE run failed, aborted.')
0565     else
0566         %
0567         Znorm = lukestructs.output.Znorm;
0568         Zcurr = lukestructs.output.Zcurr;
0569         ZP0 = lukestructs.output.ZP0;
0570         dke_out = lukestructs.output.dke_out;
0571         radialDKE = lukestructs.output.radialDKE;
0572         equilDKE = lukestructs.output.equilDKE;
0573         momentumDKE = lukestructs.output.momentumDKE;
0574         gridDKE = lukestructs.output.gridDKE;
0575         Zmomcoef = lukestructs.output.Zmomcoef;
0576         Zbouncecoef = lukestructs.output.Zbouncecoef;
0577         Zmripple = lukestructs.output.Zmripple;
0578         mksa = lukestructs.output.mksa;
0579         XXsinksource = lukestructs.output.XXsinksource;
0580         %
0581     end
0582     %
0583     clear raytracing_yp alphaphimex_jd trapz_dke_yp dmumpsmex
0584     %
0585 end
0586 %
0587 % integralle surfacique (d'une grandeur independante de theta)
0588 %  s = integrale de surface
0589 %  e = valeur a integree
0590 %  x = coordonnees normalisee
0591 %  sp = datak.equi.sp
0592 %  rhomax = datak.equi.rhomax
0593 function s=zintsurf(e,x,sp,rhomax)   
0594 
0595     s = rhomax .* trapz(x,sp .* e,2);
0596 end
cronos2luke

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
