main_dke_yp

PURPOSE ^

SYNOPSIS ^

function [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] =main_dke_yp(simul_in,dkepath,equil,dkeparam_in,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct)

DESCRIPTION ^

  Main function of the 3-D relativistic bounce-averaged electron drift kinetic equation solver 
  for the RF wave absorption problem in presence of an Ohmic electric field (electron bootstrap current consistent with non-Maxwellian distribution)
  Small banana width approximation
  Arbitrary axisymmetric plasma equilibrium (linked arbitrary tokamak magnetic equilibrium code, via mexfile)
  Plasma current, RF power absorption, runaway rate, magnetic ripple loss
  and fast electron bremsstrahlung calculations

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = ...
0002     main_dke_yp(simul_in,dkepath,equil,dkeparam_in,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in,loadstruct)
0003 %
0004 %  Main function of the 3-D relativistic bounce-averaged electron drift kinetic equation solver
0005 %  for the RF wave absorption problem in presence of an Ohmic electric field (electron bootstrap current consistent with non-Maxwellian distribution)
0006 %  Small banana width approximation
0007 %  Arbitrary axisymmetric plasma equilibrium (linked arbitrary tokamak magnetic equilibrium code, via mexfile)
0008 %  Plasma current, RF power absorption, runaway rate, magnetic ripple loss
0009 %  and fast electron bremsstrahlung calculations
0010 %
0011 %by Y.Peysson (CEA-DRFC) <yves.peysson@cea.fr> and J. Decker (CEA-DRFC) <joan.decker@cea.fr>
0012 %
0013 %
0014 if ~isstruct(simul_in),
0015     simul.id = simul_in;
0016     simul.path = './';
0017 else
0018     simul = simul_in;
0019 end
0020 %
0021 simul.timeid = timeid_jd(clock);
0022 simul.userid = userid_jd;
0023 [simul.LUKEversion,simul.matver] = LUKEversion_jd;
0024 simul.hostname = dkepath.hostname;
0025 %
0026 dke_warnings = '';
0027 %
0028 if nargin < 11,
0029     error('Wrong number of input arguments in main_dke_yp');
0030 elseif nargin == 11,
0031     loadstruct = [];
0032 end
0033 %
0034 if dkeparam_in.opt_save && ~ exist('tmp','dir'),
0035     mkdir('tmp')
0036 end
0037 %
0038 format long e;
0039 %
0040 % set temp path
0041 %
0042 mkdirmessage = 'Directory already exists.';
0043 while strcmp(mkdirmessage,'Directory already exists.'),
0044     dkepath.temppath = fullfile(dkepath.temppathroot,timeid_jd(clock,3));
0045     [dummy,mkdirmessage] = mkdir(dkepath.temppath);%Create calculation directory
0046 end
0047 %
0048 tmpfile = [simul.path,'LUKE_RESULTS_',simul.id,'_TMP.mat'];
0049 %
0050 if dkeparam_in.opt_load && exist(tmpfile,'file'),
0051     %
0052     load(tmpfile);
0053     info_dke_yp(2,['File ',tmpfile,' reloaded.']);
0054     %
0055     itn0 = itn;
0056     %
0057 else
0058     %
0059     % Load data if existing in 'loadstruct' structure or recalculate them otherwise
0060     %
0061     [dkeparam_out,radialDKE,equilDKE,momentumDKE,gridDKE,mksa,Zmomcoef,Zbouncecoef,...
0062         XXsinksource,Zmripple,XXSavalanches_norm,xnrem_init,xnrep_init,xrr_out_norm,...
0063         ZXXD_a,ZXXF_a,...
0064         ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,XXfM,XXILor,...
0065         ZXXD_e_init,ZXXF_e_init,ZXXD_e_tp,ZXXF_e_tp,xepsi_init,xEfield_validity,...
0066         ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,...
0067         waveparam,xyprop_dke,xyP0_2piRp_mod,xyP0_2piRp_mod_coll,xyP0_2piRp_mod_noabs,yb,yP0_2piRp,xys,...
0068         ZXYD_rf,ZXYF_rf,ZXYD_rf_tp,ZXYF_rf_tp,gridindex_rf,...
0069         XXfinit]...
0070          = buildlukestruct_yp(loadstruct,equil,dkeparam_in,dkedisplay,dkepath,ohm,waves,transpfaste,ripple,XXf0_in,XXsinksource_in);
0071     %
0072     time0 = clock;
0073     %
0074     if strcmp(waveparam.model,'RT') && ~isempty(waveparam.bomega_rf),%initialization of wave parameters
0075         %
0076         delta_rt = dkeparam_out.delta_rt;%smoothering parameter of quasilinear iterations
0077         %
0078         zS = raypath_jd(xys,xyprop_dke,yb,radialDKE.r_dke);
0079         nit_rf = dkeparam_out.nit_rf;
0080         if nit_rf > 1 && dkedisplay.display_mode >= 1,
0081             Fig_QLIter = figbyname_yp('Quasilinear iteration');
0082             if dkedisplay.display_mode >= 2,
0083                 Fig_Pray = figbyname_yp('Power absorption along rays');
0084             end
0085         end
0086     else
0087         nit_rf = 1;
0088         redfac = 1;
0089         xyP0_2piRp_mod = ones(gridindex_rf.nr,gridindex_rf.ny);
0090     end 
0091     %
0092     FPtime = dkeparam_out.tn;
0093     nit = length(FPtime);
0094     dtn_list = diff([0,FPtime]);
0095     %
0096     if nit > 1 && dkedisplay.display_mode >= 1,
0097         Fig_TimeIter = figure('Name','Time evolution');        
0098     end
0099     %
0100     residue_rf = cell(1,nit);
0101     normf0 = cell(1,nit);
0102     curr0 = cell(1,nit);
0103     xcurr_rem = cell(1,nit);
0104     xcurr_rep = cell(1,nit);
0105     xcurr = cell(1,nit);
0106     residu_f = cell(1,nit);
0107     XXf0 = cell(1,nit);
0108     zP_0_2piRp_mod_old = cell(nit_rf,nit);
0109     xepsi = cell(1,nit);
0110     xnorm_0 = cell(1,nit);
0111     xnorm_rem = cell(1,nit);
0112     xnorm_rep = cell(1,nit);
0113     xnorm_rim = cell(1,nit);
0114     xnorm_rip = cell(1,nit);
0115     XXRintmask = cell(1,nit);
0116     %
0117     if dkeparam_out.timevol == 1,
0118     %     ZP0 = cell(1,nit);
0119     %     Zf0_interp = cell(1,nit);
0120     %     Zcurr = cell(1,nit);
0121     %     Znorm = cell(1,nit);
0122     %     Zmom = cell(1,nit);
0123         XxRR_fsav = cell(1,nit);
0124         xMRR_flux = cell(1,nit);
0125         xMRR_tau = cell(1,nit);
0126         xMRR_power_flux = cell(1,nit);
0127         xMRR_power_tau = cell(1,nit);
0128         xRRm_fsav = cell(1,nit);
0129         xRRp_fsav = cell(1,nit);
0130         xne_norm_out = cell(1,nit);
0131         xTe_norm_out = cell(1,nit);
0132         xft_fsav = cell(1,nit);
0133         xfteff_fsav = cell(1,nit);        
0134         %
0135     elseif dkeparam_out.timevol == 2,
0136     %     ZP0 = cell(nit_rf,nit);
0137     %     Zf0_interp = cell(nit_rf,nit);
0138     %     Zcurr = cell(nit_rf,nit);
0139     %     Znorm = cell(nit_rf,nit);
0140     %     Zmom = cell(nit_rf,nit);
0141         XxRR_fsav = cell(nit_rf,nit);
0142         xMRR_flux = cell(nit_rf,nit);
0143         xMRR_tau = cell(nit_rf,nit);
0144         xMRR_power_flux = cell(nit_rf,nit);
0145         xMRR_power_tau = cell(nit_rf,nit);
0146         xRRm_fsav = cell(nit_rf,nit);
0147         xRRp_fsav = cell(nit_rf,nit);
0148         xne_norm_out = cell(nit_rf,nit);
0149         xTe_norm_out = cell(nit_rf,nit);
0150         xft_fsav = cell(nit_rf,nit);
0151         xfteff_fsav = cell(nit_rf,nit);
0152         %
0153     end
0154     %
0155     itn0 = 1;
0156     XXfin = XXfinit;
0157     xepsi_in = xepsi_init;
0158     %
0159     xqtilde = Zbouncecoef.xqtilde(gridDKE.rdke);
0160     xqhat = Zbouncecoef.xqhat(gridDKE.rdke);
0161     xq = Zbouncecoef.xq(gridDKE.rdke);
0162     xqbar = Zbouncecoef.xqbar(gridDKE.rdke);
0163     XXlambda = Zbouncecoef.XXlambda(:,:,gridDKE.rdke);
0164     %
0165     dpn = gridDKE.dpn;
0166     dmhu = gridDKE.dmhu;
0167     XXpn2 = gridDKE.XXpn2(:,:,gridDKE.rdke);
0168     XXpn = gridDKE.XXpn(:,:,gridDKE.rdke);
0169     XXmhu = gridDKE.XXmhu(:,:,gridDKE.rdke);
0170     %
0171     XXgamma = Zmomcoef.XXgamma(:,:,gridDKE.rdke);
0172     %
0173     % Internal RE domain
0174     %
0175     % TODO : - calculate the flux through the RR internal boundary
0176     %
0177     if isfield(dkeparam_out,'pnmin_S') && isscalar(dkeparam_out.pnmin_S) && ~isnan(dkeparam_out.pnmin_S),% Rfrac calculated from dkeparm.pnmin_S (see Adam's work)
0178         %
0179         % partial density above the minimum energy dkeparam.pnmin_S (typically used in runaway studies)
0180         %
0181         XXRintmask_init = gridDKE.XXpn > dkeparam_out.pnmin_S;
0182     else
0183         %
0184         % Rfrac calculated from internal RR boundary based on force balance
0185         %
0186         XXRintmask_init = ZXXF_c.p_ij + ZXXF_e_init.p_ij + ZXXF_s.p_ij > 0;% internal runaways from force balance
0187     end
0188     %
0189     % initial f0 normalized density
0190     %
0191     xnorm_0_init = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXfin(:,:,gridDKE.rdke).*XXpn2,1),1).';
0192     %
0193     % initial internal source runaway normalized density
0194     %
0195     XXRmaskmin0_KOm = gridDKE.XXpn > dkeparam_out.pnmin0_KO & gridDKE.XXmhu < 0;
0196     XXRmaskmin0_KOp = gridDKE.XXpn > dkeparam_out.pnmin0_KO & gridDKE.XXmhu > 0;
0197     %
0198     xnorm_rim_init = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXRmaskmin0_KOm(:,:,gridDKE.rdke).*XXfin(:,:,gridDKE.rdke).*XXpn2,1),1).';
0199     xnorm_rip_init = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXRmaskmin0_KOp(:,:,gridDKE.rdke).*XXfin(:,:,gridDKE.rdke).*XXpn2,1),1).';
0200     %
0201     % In frozen current mode, LUKE assumes that the current is frozen to either :
0202     %   - a prescribed scalar value (dkeparam_out.Itot_ohm) assuming a uniform electric field
0203     %   - a prescribed profile (dkeparam_out.Jtot_ohm)
0204     %   - its original value calculated from the initial distribution and number of runaways (beyond pmax)
0205     %       (It is assumed that the initial distribution function was calculated with the input value for the ohmic field)
0206     %
0207     if dkeparam_out.nit_ohm > 1 && isnan(dkeparam_out.Itot_ohm),% frozen current mode with current profile specified
0208         %
0209         if ~any(isnan(dkeparam_out.Jtot_ohm)),% prescribed current specified by dkeparam_out.Jtot_ohm on the grid dkeparam_out.rho_ohm
0210             %
0211             if length(dkeparam_out.rho_ohm) == 1,
0212                 xcurr_init = dkeparam_out.Jtot_ohm/mksa.j_ref;%Normalized current
0213             else
0214                 xcurr_init = interp1(dkeparam_out.rho_ohm,dkeparam_out.Jtot_ohm,equilDKE.xrho,'spline')/mksa.j_ref;%Normalized current
0215             end
0216             %
0217         else % prescribed current specified by xnre_init and XXfin
0218             %
0219             % runaways beyond pmax are assumed to have a parallel velocity equal to the speed of light
0220             %
0221             %
0222             % initial runaway current
0223             %
0224             xcurr_rem_init = - (xq./xqbar).*(xqhat./xqtilde).*xnrem_init/mksa.betath_ref;
0225             xcurr_rep_init =   (xq./xqbar).*(xqhat./xqtilde).*xnrep_init/mksa.betath_ref;
0226             %
0227             % initial f0 current
0228             %
0229             xcurr_0 = 2*pi*(xq./xqbar).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXfin(:,:,gridDKE.rdke).*XXpn2.*XXpn.*XXmhu./XXgamma,1),1).';
0230             %
0231             xcurr_init = xcurr_rem_init + xcurr_rep_init + xcurr_0;
0232             %
0233         end
0234     end
0235     %
0236 end
0237 %
0238 % --------------------------------- Solve the electron drift kinetic --------------------------------
0239 %
0240 % The evolution of the distribution function is calculated from time 0 to
0241 % time dkeparam.tn, by time steps dkeparam.dtn. For each time step,
0242 % nested convergence on the wave absorption and explicit collision
0243 % integral term must be obtained.
0244 %
0245 flag_exit = false;
0246 %
0247 for itn = itn0:nit,
0248     %
0249     dkeparam_out.dtn = dtn_list(itn);
0250     %
0251     % The internal and external avalanche terms are proportional to both the runaway fraction and
0252     % the bulk fraction at time itn-1
0253     %
0254     % Note : in the Rosenbluth model, we have to neglect the effect on the incident electron
0255     %
0256     [~,xne_norm,~,~,~,xbetath] = normcoef_dke_yp(mksa,equilDKE,gridDKE);
0257     %
0258     if itn == 1,   
0259         if isnan(dkeparam_out.pnmax1_KO),% count all electrons in bulk except internal runaways
0260             xnorm_b = xnorm_0_init - xnorm_rim_init;
0261         elseif dkeparam_out.pnmax1_KO == 0,% count using the value of f(0) - most logical given how the sink term is calculated
0262             xnorm_b = xne_norm.*reshape(XXfin(1,1,:)./XXfM(1,1,:),[1,gridDKE.nr]);
0263         elseif imag(dkeparam_out.pnmax1_KO) == 0,% count bulk electrons up to pnmax1_KO normalized to ref temperature
0264             xnorm_b = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*(XXpn < dkeparam_out.pnmax1_KO).*XXfin(:,:,gridDKE.rdke).*XXpn2,1),1).';
0265         else% count bulk electrons up to pnmax1_KO normalized to ref temperature
0266             xnorm_b = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*(XXpn < imag(dkeparam_out.pnmax1_KO)*repmat(reshape(xbetath,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1])/mksa.betath_ref).*XXfin(:,:,gridDKE.rdke).*XXpn2,1),1).';
0267         end
0268     else
0269         if isnan(dkeparam_out.pnmax1_KO),% count all electrons in bulk except internal runaways
0270             xnorm_b = xnorm_0{itn - 1} - xnorm_rim{itn - 1};
0271         elseif dkeparam_out.pnmax1_KO == 0,% count using the value of f(0) - most logical given how the sink term is calculated
0272             xnorm_b = xne_norm.*reshape(XXf0{itn - 1}(1,1,:)./XXfM(1,1,:),[1,gridDKE.nr]);
0273         elseif imag(dkeparam_out.pnmax1_KO) == 0,% count bulk electrons up to pnmax1_KO normalized to ref temperature
0274             xnorm_b = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*(XXpn < dkeparam_out.pnmax1_KO).*XXf0{itn-1}(:,:,gridDKE.rdke).*XXpn2,1),1).';
0275         else% count bulk electrons up to pnmax1_KO normalized to ref temperature
0276             xnorm_b = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*(XXpn < imag(dkeparam_out.pnmax1_KO)*repmat(reshape(xbetath,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1])/mksa.betath_ref).*XXf0{itn-1}(:,:,gridDKE.rdke).*XXpn2,1),1).';
0277         end
0278     end
0279     %
0280     if itn == 1,   
0281         xnorm_m = (xnrem_init + xnorm_rim_init).*xnorm_b;
0282         xnorm_p = (xnrep_init + xnorm_rip_init).*xnorm_b;
0283     else
0284         xnorm_m = (xnorm_rem{itn - 1} + xnorm_rim{itn - 1}).*xnorm_b;
0285         xnorm_p = (xnorm_rep{itn - 1} + xnorm_rip{itn - 1}).*xnorm_b;
0286     end
0287     %
0288     XXSavalanches_norm_p = XXSavalanches_norm;
0289     XXSavalanches_norm_m = flipdim(XXSavalanches_norm,2);
0290     %
0291     XXSavalanches_orig = XXSavalanches_norm_m.*repmat(reshape(xnorm_m,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1])...% multiple the avalanche source term by the runaway electron density
0292                   + XXSavalanches_norm_p.*repmat(reshape(xnorm_p,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1]);% multiple the avalanche source term by the runaway electron density
0293     %
0294     xrr_out_m = xrr_out_norm.*xnorm_m;
0295     xrr_out_p = xrr_out_norm.*xnorm_p;
0296     xrr_out = xrr_out_m + xrr_out_p;
0297     %
0298     if isnan(dkeparam_out.pnmin2_KO),% count all internal runaways
0299         if itn == 1,
0300             XXSavalanches = XXSavalanches_orig.*XXRintmask_init;
0301         else
0302             XXSavalanches = XXSavalanches_orig.*XXRintmask{itn-1};
0303         end
0304     else
0305         XXSavalanches = XXSavalanches_orig;
0306     end
0307     %
0308     if itn == 1 && dkeparam_out.finitabs == 1,%The initial distribution is used to first estimate the ray absorption
0309         conv = -1;
0310     else
0311         conv = 0;
0312     end
0313     %
0314     % Convergence on Ohmic field in frozen current mode
0315     %
0316     for it_ohm = 1:dkeparam_out.nit_ohm,
0317         %
0318         if dkeparam_out.nit_ohm > 1,
0319             if it_ohm == 1,
0320                 xepsi{itn} = zeros(1,gridDKE.nr);% first iteration without ohmic field
0321             end
0322             %
0323             ZXXD_e = structtimes_jd(ZXXD_e_init,repmat(reshape(xepsi{itn}./xepsi_init,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1]));
0324             ZXXF_e = structtimes_jd(ZXXF_e_init,repmat(reshape(xepsi{itn}./xepsi_init,[1,1,gridDKE.nr]),[gridDKE.npn,gridDKE.nmhu,1]));
0325             %ZXXD_e_tp = structtimes_jd(ZXXD_e_tp,xepsi{itn}/xepsi_init);
0326             %ZXXF_e_tp = structtimes_jd(ZXXF_e_tp,xepsi{itn}/xepsi_init);
0327             %
0328         else
0329             ZXXD_e = ZXXD_e_init;
0330             ZXXF_e = ZXXF_e_init;
0331             xepsi{itn} = xepsi_in;
0332         end
0333         %
0334         % Internal RE domain
0335         %
0336         % TODO : - calculate the flux through the RR internal boundary
0337         %
0338         if isfield(dkeparam_out,'pnmin_S') && isscalar(dkeparam_out.pnmin_S) && ~isnan(dkeparam_out.pnmin_S),% Rfrac calculated from dkeparm.pnmin_S (see Adam's work)
0339             %
0340             % partial density above the minimum energy dkeparam.pnmin_S (typically used in runaway studies)
0341             %
0342             XXRintmask{itn} = gridDKE.XXpn > dkeparam_out.pnmin_S;
0343         else
0344             %
0345             % Rfrac calculated from internal RR boundary based on force balance
0346             %
0347             XXRintmask{itn} = ZXXF_c.p_ij + ZXXF_e.p_ij + ZXXF_s.p_ij > 0;% internal runaways from force balance
0348         end
0349         %
0350         % Convergence on p=0 particle source to compensate for boundary losses
0351         %
0352         XXnormsource = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr_dke);
0353         %
0354         xRR_fsav_old = 0;
0355         %
0356         for it_norm = 1:dkeparam_out.nit_norm,
0357             %
0358             % Convergence on RF deposition with ray-tracing formalism
0359             %
0360             for it_rf = 1:nit_rf,
0361                 %
0362                 if strcmp(waveparam.model,'RT') && ~isempty(waveparam.bomega_rf),%in the case of RF waves only
0363                     %
0364                     if it_rf > 1 && conv <= 0,
0365                         xyP0_2piRp_mod = delta_rt*xyP0_2piRp_mod_old + (1 - delta_rt)*xyP0_2piRp_mod;
0366                     end
0367                     %
0368                     if dkeparam_out.redfac == 1 && nit_rf > 1,
0369                         if it_rf == 1,
0370                             redfac = max(sum(xyP0_2piRp_mod.'))/sum(max(xyP0_2piRp_mod(:,yb)));
0371                         else
0372                             redfac = redfac - 1;
0373                         end
0374                         %redfac = 10/it_rf;
0375                         %
0376                         redfac = max([redfac,1]);
0377                         %
0378                     else
0379                         redfac = 1;
0380                     end
0381                 end
0382                 %
0383                 xyP0_2piRp_mod_coll1 = xyP0_2piRp_mod_coll.*xyP0_2piRp_mod./xyP0_2piRp_mod_noabs;%adjust power remaining in fragments including NRCA and RLA
0384                 xyP0_2piRp_mod_coll1(isnan(xyP0_2piRp_mod_coll1)) = 0;
0385                 %
0386                 [ZXXD_rf,ZXXF_rf,ZXXD_rf_tp,ZXXF_rf_tp] = raysum_rf_jd(dkeparam_out,ZXYD_rf,ZXYF_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp_mod_coll1/redfac,gridindex_rf);
0387                 [Zdelta,MMXT_f_t,MMXP_f_t,MMXR_f_t,Xmask_r_edge,MMXP_tp,MMXP_g_t,MMH_tp,Zpn2_t,sm_f,sm_g,sm_tp] = matrixcalc_dke_yp(dkeparam_out,dkedisplay,gridDKE,equilDKE,Zbouncecoef,mksa,Zmripple,Zmomcoef,transpfaste,ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,XXILor,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,ZXXD_rf,ZXXF_rf,ZXXD_rf_tp,ZXXF_rf_tp,ZXXD_a,ZXXF_a);%Build the 3D DKE matrix
0388                 %
0389                 if conv >= 0,
0390                     [XXf0{itn},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,radial_mode,memoryLU_f,memoryLU_g,time_DKE,normf0{itn},curr0{itn},residu_f{itn}] = solver_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,mksa,Zpn2_t,sm_f,sm_g,sm_tp,it_rf,MMXR_f_t,MMXP_f_t,MMXT_f_t,MMXP_tp,MMXP_g_t,MMH_tp,XXILor,Xmask_r_edge,XXsinksource + XXnormsource,XXSavalanches,xrr_out,XXfM,XXfin);% Solve the DKE equation
0391                 else% case where the initial distribution is used to first estimate the ray absorption
0392                     XXf0{itn} = XXfin;
0393                     XXf0_tp = zeros(size(XXf0{itn}));
0394                     XXf0_g = zeros(size(XXf0{itn}));
0395                     XXf0_L_tp = zeros(size(XXf0{itn}));
0396                     XXf0_L_g = zeros(size(XXf0{itn}));
0397                     radial_mode = NaN;
0398                 end
0399                 %
0400                 xyP0_2piRp_mod_old = xyP0_2piRp_mod;
0401                 %
0402                 if strcmp(waveparam.model,'RT') && ~isempty(waveparam.bomega_rf),%In the case of RF waves and quasilinear iteration
0403                     if nit_rf > 1,  
0404                         %
0405                         % Calculation of the density of power absorbed for all rays
0406                         %
0407                         [ZP0_rf] = fmoments_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi{itn},waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,...
0408                             XXf0{itn},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,...
0409                             ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp_mod,gridindex_rf);%Calculate RF absorbed power of the distribution function
0410                         if dkeparam_out.dke_mode >= 1,
0411                             xyPabs_rf_dke_fsav_out = ZP0_rf.xy_rf_fsav + ZP0_rf.xy_rf_tp_fsav + ZP0_rf.xy_rf_g_fsav; % flux-surface averaged density of absorbed power.
0412                         else
0413                             xyPabs_rf_dke_fsav_out = ZP0_rf.xy_rf_fsav; % flux-surface averaged density of absorbed power.
0414                         end        
0415                         %
0416                         xyPabs_rf_dke_fsav_mksa = xyPabs_rf_dke_fsav_out*mksa.P_ref*1e6;%RF power (W)
0417                         %
0418                         % Calculation of the power flow along the rays from DKE calculations
0419                         %
0420                         [xyP0_2piRp_out,xyP0_2piRp_mod] = rfprop_dke_jd(yb,yP0_2piRp,xyprop_dke,xyP0_2piRp_mod_old,xyPabs_rf_dke_fsav_mksa,equilDKE.xdV_2piRp_dke,gridDKE.xpsin,gridDKE.xpsin_dke,gridDKE.rdke,dkeparam_out.dke_mode);
0421                         %
0422                         zP_0_2piRp_mod_old{it_rf,itn} = raypath_jd(xyP0_2piRp_mod_old,xyprop_dke,yb,radialDKE.r_dke);
0423                         zP_0_2piRp_mod = raypath_jd(xyP0_2piRp_mod,xyprop_dke,yb,radialDKE.r_dke);
0424                         zP_0_2piRp_out = raypath_jd(xyP0_2piRp_out,xyprop_dke,yb,radialDKE.r_dke);
0425                         zP_0_2piRp_mod_coll1{it_rf,itn} = raypath_jd(xyP0_2piRp_mod_coll1,xyprop_dke,yb,radialDKE.r_dke);
0426                         %
0427                         P_0_2piRp_in = sum(yP0_2piRp);
0428                         P_0_2piRp_out = 0;
0429                         %
0430                         for iray = 1:length(yb),
0431                              P_0_2piRp_out = P_0_2piRp_out + zP_0_2piRp_mod_old{it_rf,itn}{iray}(end);
0432                         end
0433                         %
0434                         abs_frac = (P_0_2piRp_in - P_0_2piRp_out)/P_0_2piRp_in;
0435                         %
0436                         for iray = 1:length(yb),
0437                             yresidue_rf_mod(iray,it_rf) = max(abs(zP_0_2piRp_mod{iray} - zP_0_2piRp_mod_old{it_rf,itn}{iray})./zP_0_2piRp_mod{iray}(1));
0438                             yresidue_rf_out(iray,it_rf) = max(abs(zP_0_2piRp_mod{iray} - zP_0_2piRp_out{iray})./zP_0_2piRp_mod{iray}(1));
0439                         end
0440                         %
0441                         %residue_rf{itn}(it_rf) = max([yresidue_rf_mod(:,it_rf);yresidue_rf_out(:,it_rf)]);
0442                         residue_rf{itn}(it_rf) = max(yresidue_rf_mod(:,it_rf));
0443                         %
0444                         % Moments of the distribution function
0445                         %
0446                         if dkeparam_out.timevol == 2,
0447                             [ZP0(it_rf,itn),Zcurr(it_rf,itn),Znorm(it_rf,itn),ZXXS,...(itn)
0448                                 XxRR_fsav{it_rf,itn},xMRR_flux{it_rf,itn},xMRR_tau{it_rf,itn},xMRR_power_flux{it_rf,itn},xMRR_power_tau{it_rf,itn},xRRm_fsav{it_rf,itn},xRRp_fsav{it_rf,itn},...
0449                                 xne_norm_out{it_rf,itn},xTe_norm_out{it_rf,itn},xft_fsav{it_rf,itn},xfteff_fsav{it_rf,itn},Zmom(it_rf,itn)] = ...
0450                                 fmoments_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi{itn},waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,...
0451                                 XXf0{itn},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,...
0452                                 ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp_mod_old,gridindex_rf,...
0453                                 ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,XXRintmask{itn},XXILor); %Calculate moments of the distribution function
0454                         end
0455                         %
0456                         if dkedisplay.display_mode >= 1,
0457                             figure(Fig_QLIter),
0458                             [ax] = graph1D_jd(1:it_rf,residue_rf{itn},0,1,'Iteration number','\DeltaP/P','Quasilinear convergence',NaN,[0,nit_rf],[dkeparam_out.prec_rf,1],'none','+','r',2);
0459                             grid  on
0460                             drawnow;                
0461             %
0462                             if dkedisplay.display_mode >= 2,
0463                                 for iray = 1:length(zS),
0464                                     figure(Fig_Pray + iray - 1);clf
0465                                     [ax] = graph1D_jd(zS{iray},zP_0_2piRp_mod_old{it_rf,itn}{iray},0,0,'','P_{LH}/2\piR_p (W/m)',['QL convergence, ray #',int2str(iray)],NaN,'','','--','none','r',2);
0466                                     [ax] = graph1D_jd(zS{iray},zP_0_2piRp_mod{iray},0,0,'','P_{LH}/2\piR_p (W/m)',['QL convergence, ray #',int2str(iray)],NaN,'','','-','none','r',2);
0467                                     legend(['Iteration #',int2str(it_rf-1)],['Iteration #',int2str(it_rf)])
0468                                     title(['Ray #',num2str(iray),', residue: ',num2str(max([yresidue_rf_mod(iray,it_rf);yresidue_rf_out(iray,it_rf)]))])
0469                                     grid on
0470                                     drawnow;
0471                                 end
0472                             end
0473                         end
0474             %
0475                         if dkedisplay.display_mode >=1,
0476                             info_dke_yp(2,'---------------------------------------------------------');
0477                             info_dke_yp(2,['Quasilinear iteration : ',int2str(it_rf),', Residu = ',num2str(residue_rf{itn}(it_rf)),', delta_rt = ',num2str(delta_rt)]);
0478                             info_dke_yp(2,'---------------------------------------------------------');
0479                         end
0480             %
0481                         if (residue_rf{itn}(it_rf) <= dkeparam_out.prec_rf && redfac == 1 && it_rf > 1 && conv == 1) || it_rf == nit_rf,
0482                             %
0483                             dke_warnings = '';               
0484                             if isfield(dkeparam_out,'abs_lim') && abs_frac < dkeparam_out.abs_lim,
0485                                 dke_warnings = 'The RF power is not fully absorbed. To bypass this error set dkeparam.abs_lim between 0 and 1';
0486                             end
0487                             %
0488                             if residue_rf{itn}(it_rf) <= dkeparam_out.prec_rf && redfac == 1 && it_rf > 1 && conv == 1,
0489                                 if dkedisplay.display_mode >=1,info_dke_yp(2,['Successfull convergence achieved for ray-tracing coupling after ',int2str(it_rf),' iterations !']);end
0490                                 break
0491                             else
0492                                 info_dke_yp(2,['WARNING: convergence for ray-tracing coupling was not achieved after ',int2str(it_rf),' iterations !']);                  
0493                             end
0494                         else       
0495                             if dkedisplay.display_mode >=1,info_dke_yp(2,['Iteration ',int2str(it_rf),', residue = ',num2str(residue_rf{itn}(it_rf))]);end       
0496                         end
0497                         %
0498                         if (residue_rf{itn}(it_rf) <= dkeparam_out.prec_rf && redfac == 1 && it_rf > 1 && conv == 0),
0499                             if dkeparam_out.postest_rf && min(min(min(ZP0_rf.xyn_rf_fsav))) < 0,
0500                                  info_dke_yp(2,['WARNING: residue below threshold but iterations continue while negative absorbed power occurs somewhere!']);                  
0501                             else
0502                                 conv = 1;
0503                             end
0504                         end        
0505                         %
0506                         if (residue_rf{itn}(it_rf) <= dkeparam_out.prec_rf && redfac == 1 && it_rf > 1 && conv == -1),
0507                             if dkedisplay.display_mode >=1,info_dke_yp(2,['Successfull convergence achieved for ray-tracing coupling on initial distribution after ',int2str(it_rf),' iterations !']);end
0508                             conv = 0;
0509                         end    
0510                         %
0511                         if 0,%dkeparam_out.finsert == 1,
0512                             XXfin = XXf0{itn};
0513                             abs_frac_rf(it_rf) = abs_frac;
0514                         end
0515                     end     
0516                 end
0517             end
0518             %
0519             % Moments of the distribution function
0520             %
0521             if dkeparam_out.timevol == 1,
0522                 [ZP0(itn),Zcurr(itn),Znorm(itn),ZXXS,...(itn)
0523                     XxRR_fsav{itn},xMRR_flux{itn},xMRR_tau{itn},xMRR_power_flux{itn},xMRR_power_tau{itn},xRRm_fsav{itn},xRRp_fsav{itn},...
0524                     xne_norm_out{itn},xTe_norm_out{itn},xft_fsav{itn},xfteff_fsav{itn},Zmom(itn)] = ...
0525                     fmoments_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi{itn},waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,...
0526                     XXf0{itn},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,...
0527                     ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp_mod_old,gridindex_rf,...
0528                     ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,XXRintmask{itn},XXILor); %Calculate moments of the distribution function
0529                 %
0530             end
0531             %
0532             % compensation term for boundary losses
0533             %
0534             if dkeparam_out.nit_norm > 1,
0535                 %
0536                 xRR_fsav = xRRm_fsav{itn} + xRRp_fsav{itn};
0537                 %
0538                 if all(abs(xRR_fsav) <= abs(dkeparam_out.thres_norm)),% negligible boundary losses
0539                     if dkedisplay.display_mode >= 1, 
0540                         disp(['Boundary losses negligible, no p=0 source term is needed at time step ',num2str(itn)])
0541                     end
0542                     break
0543                 elseif max(abs(xRR_fsav - xRR_fsav_old)./abs(xRR_fsav)) <= dkeparam_out.tol_norm,
0544                     if dkedisplay.display_mode >= 1, 
0545                         disp(['Convergence on p=0 source term reached after ',num2str(it_norm),' iterations at time step ',num2str(itn)])
0546                     end
0547                     break
0548                 end
0549                 %
0550                 XXnormsource(1,:,:) = repmat(reshape(xRR_fsav.*xqhat./xqtilde,[1,1,gridDKE.nr_dke]),[1,gridDKE.nmhu,1])./XXlambda(1,:,:)/(4*pi*gridDKE.dpn(1)*gridDKE.pn(1).^2);
0551                 %
0552                 xRR_fsav_old = xRR_fsav;
0553                 %
0554             end
0555         end
0556         %
0557         if dkeparam_out.nit_norm > 1 && any(abs(xRR_fsav) > abs(dkeparam_out.thres_norm)) && max(abs(xRR_fsav - xRR_fsav_old)./abs(xRR_fsav)) > dkeparam_out.tol_norm,
0558             disp(['WARNING : Convergence on p=0 source term not reached after ',num2str(it_norm),' iterations at time step ',num2str(itn)])
0559         end
0560         %
0561         if dkeparam_out.timevol == 1,
0562             %
0563             % normalization to compensate for numerical errors
0564             %
0565             % if imag(dkeparam.norm_mode_f) == 1, particle conservation is enforced but
0566             % boundary losses are accounted for (including secondary electrons born
0567             % outside LUKE domain, which are compensated for in solver_dke_yp).
0568             %
0569             % with real(dkeparam_out.norm_mode_f) == 0 , this should only compensate for numerical errors.
0570             % with real(dkeparam_out.norm_mode_f) == 2 , this should also compensate for boundary losses, which are canceled out in solver_dke_yp.
0571             % with real(dkeparam_out.norm_mode_f) == 1 , this option should not be used as the norm of f is enforced to the Maxwellian value
0572             %
0573             if Znorm(itn).x_0_fsav < 0,
0574                 disp('WARNING : negative norm of the electron distribution')
0575                 flag_exit = true;
0576             end
0577             %
0578             if ~flag_exit && imag(dkeparam_out.norm_mode_f) == 1,
0579                 %
0580                 if itn == 1,
0581                     xnorm_0_last = xnorm_0_init;% LUKE edf normalization at previous step
0582                 else
0583                     xnorm_0_last = xnorm_0{itn - 1};% LUKE edf normalization at previous step
0584                 end
0585                 %
0586                 xnorm_0_corr = xnorm_0_last - (xrr_out_m + xrr_out_p)*dtn_list(itn) + xne_norm_out{itn};% subtract avalanches RE created outside LUKE domain, plus other sources/sinks
0587                 dxnorm_flux = (xRRm_fsav{itn} + xRRp_fsav{itn} - xRR_fsav_old)*dtn_list(itn);% RE leaving LUKE domain (possibly corrected by compensation)
0588                 %
0589                 % this expressions derive from the fact that xRR_fsav{itn} is itself proportional to f
0590                 %
0591                 % it brings first order corrections to the approximate expression :
0592                 %   Znorm(itn).xnorm_fac = Znorm(itn).x_0_fsav./(xnorm_0_corr - dxnorm_flux);
0593                 %
0594 %                 xnorm_fac = Znorm(itn).x_0_fsav./xnorm_0_corr.*(1 + Znorm(itn).x_0_fsav.*dxnorm_flux./xnorm_0_corr.^2);
0595                 xnorm_fac = (Znorm(itn).x_0_fsav + dxnorm_flux)./xnorm_0_corr;
0596                 %
0597                 % renormalize f moments accordingly
0598                 %
0599                 [XXf0{itn},ZP0(itn),Zcurr(itn),Znorm(itn),ZXXS,...
0600                         XxRR_fsav{itn},xMRR_flux{itn},xMRR_tau{itn},xMRR_power_flux{itn},xMRR_power_tau{itn},xRRm_fsav{itn},xRRp_fsav{itn},xTe_norm_out{itn},Zmom(itn)] = ...
0601                     renorm_fmoments_jd(xnorm_fac,equilDKE,mksa,...
0602                         XXf0{itn},ZP0(itn),Zcurr(itn),Znorm(itn),ZXXS,...
0603                         XxRR_fsav{itn},xMRR_flux{itn},xMRR_tau{itn},xMRR_power_flux{itn},xMRR_power_tau{itn},xRRm_fsav{itn},xRRp_fsav{itn},xTe_norm_out{itn},Zmom(itn));
0604                 %
0605             end
0606             %
0607             xnorm_0{itn} = Znorm(itn).x_0_fsav;
0608             %
0609             % the external runaway increment is the sum of :
0610             %   - the runaways created outside pmax
0611             %   - the runaways passing through pmax
0612             %
0613             dxnorm_rem = (xrr_out_m + xRRm_fsav{itn})*dtn_list(itn);
0614             dxnorm_rep = (xrr_out_p + xRRp_fsav{itn})*dtn_list(itn);
0615             %
0616 %             xnorm_ri{itn} = Znorm(itn).x_Rfrac_0_fsav;
0617             xnorm_rim{itn} = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXRmaskmin0_KOm(:,:,gridDKE.rdke).*XXf0{itn}(:,:,gridDKE.rdke).*XXpn2,1),1).';
0618             xnorm_rip{itn} = 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXRmaskmin0_KOp(:,:,gridDKE.rdke).*XXf0{itn}(:,:,gridDKE.rdke).*XXpn2,1),1).';
0619             %
0620         else
0621             %
0622             dxnorm_rem = zeros(1,gridDKE.nr);
0623             dxnorm_rep = zeros(1,gridDKE.nr);
0624             %
0625 %             xnorm_ri{itn} = 0;
0626             xnorm_rim{itn} = 0;
0627             xnorm_rip{itn} = 0;
0628             %
0629         end
0630         %
0631         if itn == 1,
0632             xnorm_rem{itn} = xnrem_init + dxnorm_rem;
0633             xnorm_rep{itn} = xnrep_init + dxnorm_rep;
0634         else
0635             xnorm_rem{itn} = xnorm_rem{itn - 1} + dxnorm_rem;
0636             xnorm_rep{itn} = xnorm_rep{itn - 1} + dxnorm_rep;
0637         end                
0638         %
0639         % runaway current
0640         %
0641         xcurr_rem{itn} = - (xq./xqbar).*(xqhat./xqtilde).*xnorm_rem{itn}/mksa.betath_ref;
0642         xcurr_rep{itn} =   (xq./xqbar).*(xqhat./xqtilde).*xnorm_rep{itn}/mksa.betath_ref;
0643         %
0644         % initial f0 current
0645         %
0646         xcurr_0 = 2*pi*(xq./xqbar).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0{itn}(:,:,gridDKE.rdke).*XXpn2.*XXpn.*XXmhu./XXgamma,1),1).';
0647         %
0648         xcurr{itn} = xcurr_rem{itn} + xcurr_rep{itn} + xcurr_0;% total current
0649         %
0650         % adjust ohmic field in frozen current mode
0651         %
0652         if dkeparam_out.nit_ohm > 1,% frozen current mode
0653             %
0654             if isnan(dkeparam_out.Itot_ohm),% convergence on the current profile
0655                 %
0656                 if imag(dkeparam_out.tol_ohm) > 0,% absolute threshold
0657                     tol_ohm = imag(dkeparam_out.tol_ohm);
0658                 else%relative threshold
0659                     tol_ohm = dkeparam_out.tol_ohm*max(abs(xcurr_init));
0660                 end
0661                 %
0662                 if max(abs(xcurr{itn} - xcurr_init)) <= tol_ohm,% *dkeparam_out.dtn,
0663                     if dkedisplay.display_mode >= 1,
0664                         disp(['Convergence reached after ',num2str(it_ohm),' iterations at time step ',num2str(itn)])
0665                     end
0666                     break
0667                 elseif it_ohm == 1,% first iteration without electric field
0668                     %
0669                     xcurr_in = xcurr{itn};
0670                     xepsi{itn} = xepsi_in;% previous converged field is first guess for new time
0671                     %
0672                 elseif it_ohm == dkeparam_out.nit_ohm,
0673                     disp(['WARNING : Convergence not reached after ',num2str(it_ohm),' iterations at time step ',num2str(itn)])
0674                 else
0675                     %
0676                     xepsi{itn} = xepsi{itn}.*(1 - (xcurr{itn} - xcurr_init)./(xcurr{itn} - xcurr_in));
0677                     %
0678                     if gridDKE.nr > 1,
0679                         xepsi{itn}(end) = xepsi{itn}(end - 1);% to avoid edge pb
0680                     end
0681                     %
0682                 end
0683                 %
0684             else% convergence on the total current
0685                 %
0686                 I_tot = sum(xcurr{itn}.*equilDKE.xdA_dke*mksa.j_ref);% total driven current
0687                 %
0688                 if it_ohm == 1,% first iteration without electric field
0689                     %
0690                     I_tot_in = I_tot;
0691                     xepsi{itn} = xepsi_in;% previous converged field is first guess for new time
0692                     %
0693                 elseif abs(I_tot - dkeparam_out.Itot_ohm)/abs(dkeparam_out.Itot_ohm) <= dkeparam_out.tol_ohm,% *dkeparam_out.dtn,
0694                     if dkedisplay.display_mode >= 1,
0695                         disp(['Convergence reached after ',num2str(it_ohm),' iterations at time step ',num2str(itn)])
0696                     end
0697                     break
0698                 elseif it_ohm == dkeparam_out.nit_ohm,
0699                     disp(['WARNING : Convergence not reached after ',num2str(it_ohm),' iterations at time step ',num2str(itn)])
0700                 else
0701                     %
0702                     xepsi{itn} = xepsi{itn}*(1 - (I_tot - dkeparam_out.Itot_ohm)./(I_tot - I_tot_in));
0703                     %
0704                 end                
0705                 %
0706             end
0707             %
0708         end
0709     end
0710     %
0711     xepsi_in = xepsi{itn};
0712     %
0713     if nit > 1,
0714         if dkedisplay.display_mode >= 1,
0715             figure(Fig_TimeIter),
0716             if length(equilDKE.xrho) == 1,
0717                 graph1D_jd(FPtime(itn),curr0{itn}(end),0,0,'FP time (\tau_c)','J (norm.)','FP Evolution',NaN,[0,max(FPtime)],NaN,'none','+','r',2);
0718             elseif dkeparam_out.timevol == 1,
0719                 graph1D_jd(FPtime(itn),Zcurr(itn).I_tot,0,0,'FP time (\tau_c)','I (MA)','FP Evolution',NaN,[0,max(FPtime)],NaN,'none','+','r',2);
0720             elseif dkeparam_out.timevol == 2,
0721                 graph1D_jd(FPtime(itn),Zcurr(it_rf,itn).I_tot,0,0,'FP time (\tau_c)','I (MA)','FP Evolution',NaN,[0,max(FPtime)],NaN,'none','+','r',2);
0722             else
0723                 graph1D_jd(FPtime(itn),max(abs(curr0{itn}(:,end))),0,0,'FP time (\tau_c)','max(|J|) (norm.)','FP Evolution',NaN,[0,max(FPtime)],NaN,'none','+','r',2);
0724             end
0725             grid on
0726             drawnow;                
0727         else
0728             disp(['time step ',num2str(itn),'/',num2str(nit),' completed.'])
0729         end
0730     end
0731     %
0732     if dkeparam_in.opt_save,
0733         save(tmpfile);
0734     end
0735     %
0736     XXfin = XXf0{itn};
0737     %
0738     if flag_exit,
0739         if isfield(dkeparam_out,'keyboard') && dkeparam_out.keyboard,
0740             keyboard;
0741         else
0742             break;
0743         end
0744     end
0745     %
0746 end
0747 %
0748 if nit > 1 && dkedisplay.display_mode >= 1,
0749     close(Fig_TimeIter)
0750 end
0751 %
0752 % Moments of the distribution function
0753 %
0754 if dkeparam_out.timevol == 0,
0755     %
0756     XXf0 = XXf0(nit);
0757     residue_rf = residue_rf(nit);
0758     xepsi = xepsi(nit);
0759     %
0760     [ZP0,Zcurr,Znorm,ZXXS,...
0761         XxRR_fsav,xMRR_flux,xMRR_tau,xMRR_power_flux,xMRR_power_tau,xRRm_fsav,xRRp_fsav,...
0762         xne_norm_out,xTe_norm_out,xft_fsav,xfteff_fsav,Zmom] = ...
0763         fmoments_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,xepsi{1},waveparam,radial_mode,Zdelta,XXsinksource,XXSavalanches,...
0764         XXf0{1},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,XXfM,...
0765         ZXYD_rf,ZXYD_rf_tp,ZXYF_rf_tp,xyP0_2piRp_mod_old,gridindex_rf,...
0766         ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp,XXRintmask{1},XXILor);%Calculate moments of the distribution function
0767     %
0768 end
0769 %   RF wave quasilinear diffusion
0770 %
0771 if ~isstruct(waveparam) && ~isempty(waveparam.bomega_rf),%In the case of RF waves and quasilinear iteration
0772     [qe,me] = pc_dke_yp;%Universal physics constants
0773     blist_lh = find(round(waveparam(1,:,1)./(qe*min(equilDKE.xB0)/me)) == 0);
0774     nb_lh = length(blist_lh);
0775     waveparam_lh = waveparam(:,blist_lh,:);
0776     blist_ec = find(round(waveparam(1,:,1)./(qe*min(equilDKE.xB0)/me)) > 0);
0777     nb_ec = length(blist_ec);
0778     waveparam_ec = waveparam(:,blist_ec,:);
0779 %
0780 %For displaying RF quasilinear domain only
0781 %
0782    if nb_ec > 0,
0783       xDperp_cy_ref = waveparam_ec(:,1,4);
0784    end
0785    if nb_lh > 0,
0786       xD0_lh_ref = waveparam_lh(:,1,4);
0787       xvparmin_lh = 1./waveparam_lh(:,1,6)/mksa.betath_ref;
0788       xvparmax_lh = 1./waveparam_lh(:,1,5)/mksa.betath_ref;
0789    end
0790 else 
0791     nb_lh = 0;
0792     nb_ec = 0;
0793 end
0794 %
0795 time_out = etime(clock,time0);
0796 memoryLU_f_out = memoryLU_f;
0797 memoryLU_g_out = memoryLU_g;
0798 %
0799 if dkedisplay.display_mode >= 1,
0800     if dkeparam_out.dke_mode == 0,
0801         info_dke_yp(2,['Full elapsed time for solving the electron Fokker-Planck equation: ',num2str(time_out),' (s)']) 
0802     else
0803         info_dke_yp(2,['Full elapsed time for solving the electron Drift Kinetic equation: ',num2str(time_out),' (s)']) 
0804     end
0805     disp('-------------------------------------------------------------------------------------------------------------------');
0806 end
0807 %
0808 %   Maxwellian bootstrap current calculations (for synergy calculations), no RF waves
0809 %
0810 if dkeparam_out.syn_mode == 1,
0811     time0_bc = clock;
0812     %
0813     waveparam_bc = [];
0814     [ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,xepsi,xEfield_validity] = efield_dke_jd(dkepath,dkeparam_out,dkedisplay,equilDKE,mksa,gridDKE,Zmomcoef,Zbouncecoef,'',ZXXF_c);%No Ohmic electric field
0815     %
0816     [ZYXXD_rf,ZYXXF_rf,ZYXXD_rf_tp,ZYXXF_rf_tp,ylist_rf] = rfdiff_dke_jd(dkeparam_out,dkedisplay,waveparam_bc,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,mksa);% RF wave flux coefficients
0817     [MMXP_f_t,MMXT_f_t,MMXP_tp,MMXP_g_t,MMH_tp,Zpn2_t,sm_f,sm_g,sm_tp] = matrixcalc_dke_yp(dkeparam_out,dkedisplay,gridDKE,equilDKE,Zbouncecoef,mksa,Zmripple,ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,XXILor,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZYXXD_rf,ZYXXF_rf,ZYXXD_rf_tp,ZYXXF_rf_tp);%Build the 3D DKE matrix
0818     [XXf0_bc,XXf0_bc_tp,XXf0_bc_g,XXf0_bc_L_tp,XXf0_bc_L_g,memoryLU_bc_f,memoryLU_bc_g] = solver_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,mksa,Zpn2_t,sm_f,sm_g,sm_tp,it_rf,MMXR_f_t,MMXP_f_t,MMXT_f_t,MMXP_tp,MMXP_g_t,MMH_tp,XXILor,XXsinksource,XXSavalanches,XXfinit);% Solve the DKE equation
0819     [ZP0_bc,Zcurr_bc,Znorm_bc,ZXXS_bc,ZYXXD_bc_rf,XxRR_bc_fsav,xMRR_flux_bc,xMRR_tau_bc,xMRR_power_flux_bc,xMRR_power_tau_bc,xRRm_bc_fsav,xRRp_bc_fsav,xne_norm_bc_out_bc,xTe_norm_bc_out,xft_bc_fsav,xfteff_bc_fsav] = fmoments_dke_yp(dkepath,dkeparam_out,dkedisplay,equilDKE,gridDKE,Zmomcoef,Zbouncecoef,mksa,[],radial_mode,Zdelta,XXSavalanches,XXf0{nit},XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,ZXXD_rf,ZXXD_rf_tp,ZXXF_rf_tp,ZXXD_c,ZXXF_c,ZXXD_c_tp,ZXXF_c_tp,ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,ZXXD_s,ZXXF_s,ZXXD_s_tp,ZXXF_s_tp);%Calculate moments of the distribution function
0820 %
0821     time_bc_out = etime(clock,time0_bc);
0822     memoryLU_bc_f_out = memoryLU_bc_f;
0823     memoryLU_bc_g_out = memoryLU_bc_g;
0824 %
0825     if dkedisplay.display_mode >= 1,
0826         if dke_mode == 0,
0827             info_dke_yp(2,['Full elapsed time for solving the electron Fokker-Planck equation (Bootstrap Calc. Only): ',num2str(time_bc_out),' (s)']) 
0828         else
0829             info_dke_yp(2,['Full elapsed time for solving the electron Drift Kinetic equation (Bootstrap Calc. Only): ',num2str(time_bc_out),' (s)']) 
0830         end
0831         disp('-------------------------------------------------------------------------------------------------------------------');
0832     end
0833 %
0834 end
0835 %
0836 % DKE data output structure
0837 %
0838 dke_out.XXfinit = XXfinit;
0839 dke_out.XXf0 = XXf0;
0840 dke_out.XXfM = XXfM;
0841 dke_out.FPtime = FPtime;
0842 dke_out.normf0 = normf0;
0843 dke_out.curr0 = curr0;
0844 dke_out.xcurr = xcurr;
0845 dke_out.xcurr_rem = xcurr_rem;
0846 dke_out.xcurr_rep = xcurr_rep;
0847 dke_out.residu_f = residu_f;
0848 dke_out.ZXXS = ZXXS;
0849 %
0850 dke_out.ZXXD_c = ZXXD_c;
0851 dke_out.ZXXF_c = ZXXF_c;
0852 dke_out.XXILor = XXILor;
0853 %
0854 dke_out.ZXXD_e = ZXXD_e;
0855 dke_out.ZXXF_e = ZXXF_e;
0856 dke_out.xepsi_init = xepsi_init;
0857 dke_out.xepsi = xepsi;
0858 dke_out.xEfield_validity = xEfield_validity;
0859 %
0860 dke_out.waveparam = waveparam;
0861 dke_out.xyprop_dke = xyprop_dke;
0862 dke_out.xyP0_2piRp_mod = xyP0_2piRp_mod_old;
0863 if exist('xyP0_2piRp_mod_coll1')
0864     dke_out.xyP0_2piRp_mod_coll = xyP0_2piRp_mod_coll1;
0865 end
0866 dke_out.yb = yb;
0867 dke_out.yP0_2piRp = yP0_2piRp;
0868 dke_out.xys = xys;
0869 %
0870 dke_out.ZXXD_s = ZXXD_s;
0871 dke_out.ZXXF_s = ZXXF_s;
0872 %
0873 dke_out.ZXXD_rf = ZXXD_rf;
0874 dke_out.ZXXF_rf = ZXXF_rf;
0875 %
0876 if dkeparam_out.opt_struct > 1,% save heavy structures containing individual ray element contributions to Dql
0877     dke_out.ZXYD_rf = ZXYD_rf;
0878     dke_out.ZXYF_rf = ZXYF_rf;
0879     dke_out.gridindex_rf = gridindex_rf;
0880 end
0881 %
0882 dke_out.ZXXD_a = ZXXD_a;
0883 dke_out.ZXXF_a = ZXXF_a;
0884 %
0885 dke_out.XxRR_fsav = XxRR_fsav;
0886 dke_out.xMRR_flux = xMRR_flux;
0887 dke_out.xMRR_tau = xMRR_tau;
0888 dke_out.xMRR_power_flux = xMRR_power_flux;
0889 dke_out.xMRR_power_tau = xMRR_power_tau;
0890 dke_out.xTe_norm_out = xTe_norm_out;
0891 dke_out.xne_norm_out = xne_norm_out;
0892 dke_out.xft_fsav = xft_fsav;
0893 dke_out.xfteff_fsav = xfteff_fsav;
0894 dke_out.XXRintmask = XXRintmask;
0895 dke_out.Zmom = Zmom;
0896 dke_out.memoryLU_f_out = memoryLU_f_out;
0897 dke_out.memoryLU_g_out = memoryLU_g_out;
0898 dke_out.time_out = time_out;
0899 dke_out.time_DKE = time_DKE;
0900 %
0901 if strcmp(waveparam.model,'RT') && ~isempty(waveparam.bomega_rf) && nit_rf > 1,
0902     dke_out.yresidue_rf_mod = yresidue_rf_mod;
0903     dke_out.yresidue_rf_out = yresidue_rf_out;
0904     dke_out.residue_rf = residue_rf;
0905     dke_out.zS = zS;
0906     if dkeparam_in.timevol == 0,
0907         dke_out.zP_0_2piRp_mod = zP_0_2piRp_mod_old(:,end);%last time only
0908     else
0909         dke_out.zP_0_2piRp_mod = zP_0_2piRp_mod_old;
0910     end
0911     if exist('zP_0_2piRp_mod_coll1'),
0912         if dkeparam_in.timevol == 0,
0913             dke_out.zP_0_2piRp_mod_coll = zP_0_2piRp_mod_coll1(:,end);%last time only
0914         else   
0915             dke_out.zP_0_2piRp_mod_coll = zP_0_2piRp_mod_coll1;
0916         end
0917     end
0918     %
0919     %waves = alphaphi_luke_jd(waves,dke_out,equilDKE);% calculates the QL ray absorption
0920     %
0921 end
0922 %
0923 dke_out.xRRm_fsav = xRRm_fsav;
0924 dke_out.xRRp_fsav = xRRp_fsav;
0925 %
0926 dke_out.xnorm_rem = xnorm_rem;
0927 dke_out.xnorm_rep = xnorm_rep;
0928 dke_out.xnrem_init = xnrem_init;
0929 dke_out.xnrep_init = xnrep_init;
0930 dke_out.xrr_out_norm = xrr_out_norm;
0931 dke_out.Zmripple = Zmripple;
0932 %
0933 dke_out.equil = equil;
0934 dke_out.waves = waves;
0935 dke_out.simul = simul;
0936 %
0937 if exist('abs_frac_rf','var'),
0938     dke_out.abs_frac_rf = abs_frac_rf;
0939 end
0940 %
0941 if dkeparam_out.syn_mode == 1,
0942     dke_out.XXf0_bc = XXf0_bc;
0943     dke_out.XXf0_bc_tp = XXf0_bc_tp;
0944     dke_out.XXf0_bc_g = XXf0_bc_g;
0945     dke_out.XXf0_bc_L_tp = XXf0_bc_L_tp;
0946     dke_out.XXf0_bc_L_g = XXf0_bc_L_g;
0947     dke_out.Znorm_bc = Znorm_bc;
0948     dke_out.Zcurr_bc = Zcurr_bc;
0949     dke_out.ZP0_bc = ZP0_bc;
0950     dke_out.ZXXS_bc = ZXXS_bc;
0951     dke_out.ZYXXD_bc_rf = ZYXXD_bc_rf;
0952     dke_out.XxRR_bc_fsav = XxRR_bc_fsav;
0953     dke_out.xMRR_flux_bc = xMRR_flux_bc;
0954     dke_out.xMRR_tau_bc = xMRR_tau_bc;
0955     dke_out.xMRR_power_flux_bc = xMRR_power_flux_bc;
0956     dke_out.xMRR_power_tau_bc = xMRR_power_tau_bc;
0957     dke_out.xne_norm_bc_out = xne_norm_bc_out;
0958     dke_out.xTe_norm_bc_out = xTe_norm_bc_out;
0959     dke_out.xft_bc_fsav = xft_bc_fsav;
0960     dke_out.xfteff_bc_fsav = xfteff_bc_fsav;
0961     dke_out.memoryLU_bc_g_out = memoryLU_bc_g_out;
0962     dke_out.memoryLU_bc_f_out = memoryLU_bc_f_out;
0963     dke_out.time_bc_out = time_bc_out;
0964     dke_out.xRRm_bc_fsav = xRRm_bc_fsav;
0965     dke_out.xRRp_bc_fsav = xRRp_bc_fsav;
0966 else
0967     dke_out.XXf0_bc = [];
0968     dke_out.XXf0_bc_tp = [];
0969     dke_out.XXf0_bc_g = [];
0970     dke_out.XXf0_bc_L_tp = [];
0971     dke_out.XXf0_bc_L_g = [];
0972     dke_out.Znorm_bc = [];
0973     dke_out.Zcurr_bc = [];
0974     dke_out.ZP0_bc = [];
0975     dke_out.ZXXS_bc = [];
0976     dke_out.ZYXXD_bc_rf = [];
0977     dke_out.XxRR_bc_fsav = [];
0978     dke_out.xMRR_flux_bc = [];
0979     dke_out.xMRR_tau_bc = [];
0980     dke_out.xMRR_power_flux_bc = [];
0981     dke_out.xMRR_power_tau_bc = [];    
0982     dke_out.xne_norm_bc_out = [];
0983     dke_out.xTe_norm_bc_out = [];
0984     dke_out.xft_bc_fsav = [];
0985     dke_out.xfteff_bc_fsav = [];
0986     dke_out.memoryLU_bc_g_out = [];
0987     dke_out.memoryLU_bc_f_out = [];
0988     dke_out.time_bc_out = [];
0989     dke_out.xRRm_bc_fsav = [];
0990     dke_out.xRRp_bc_fsav = [];
0991 end
0992 %
0993 if exist(dkepath.temppath,'dir')
0994     try
0995        rmdir(dkepath.temppath);%Remove calculation directory
0996     catch
0997        disp(['WARNING: The calculation directory:',dkepath.temppath,' has NOT been removed.']),
0998     end
0999 end
1000 %
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.