solver_dke_yp

   Self-consistent 3-D relativistic solver of the drift kinetic equation in momentum space with Beliaev-Budker collision integrals 
   Includes toroidal magnetic losses for the Fokker-Planck equation only for circular plasma cross-section (TS)
   Fully implicit or Crank-Nicholson time differencing schemes
   LH and EC wave quasilinear diffusion calculations + ohmic electric field 
   IMPORTANT: All input quantities are defined on the distribution function grid: half grid

   Version 4: non-uniform momentum and pitch-angle grids + fully implicit 3D calculations

   See references: 

       - C. Karney, Comput. Phys. Rep. 4 (1986) 183.
       - B. Braams, C. Karney, Phys.fluids B1 (1989) 1355
       - M. Shoucri and I. Shkarofsky, Comp. Phys. Comm. 82 (1994) 287
       - I. Shkarofsky and M. Shoucri, Nuclear Fusion, 37 (1997) 539
       - Y. Peysson and M. Shoucri, Comp. Phys. Comm. 109 (1998) 55
       - S.D. Schultz, A. Bers and A.B. Rahm, AIP Proc. of the 12th Conf. Radio Frequency Power in Plasmas, Savannah (USA), 403 1997,327
       - M. Ju et al., Phys. Plasmas, 9 (2002) 493

   (see test_dke_1yp.m)

   Input:

       - betath: Normalized reference thermal velocity (vth/c) [1,1]
       - Rp_norm: Normalized plasma major radius (ap) [1,1]
       - xTe_norm: normalized electron temperature (Te/Te_ref) [1,nr]
       - xne_norm: normalized electron density (ne/ne_ref) [1,nr]
       - Zi: ion charges [1,m]
       - mi: ion masses (u.m.a.) [1,m]
       - xni_norm: normalized ion density (ni/ne_ref) [m,nr]
       - pn: normalized electron momentum space grid (pth_ref) [1,np]
       - mhu: cos(theta) angular grid [1,nmhu]
       - XXDlh: LH normalized quasi-linear diffusion coefficient (nhuth_ref*pth_ref^2) [np,nmhu,nr]
       - XXDcy: EC normalized quasi-linear diffusion coefficient (nhuth_ref*pth_ref^2) [np,nmhu,nr]
       - param_cy: *******************normalized EC frequency [1,nr]
       - XXsinksource: normalized source/loss frequency (nhuth_ref) [np,nmhu,nr]
       - xepsi_fsav: flux surface averaged normalized ohmic electric field (pth_ref*nhuth_ref/e) [1,nr]
       - rdke: list of radial positions indexes at which FP or DKE are solved [1,nr_dke].  
       - bounce: First argument: bounce averaged option, (0) no trapping, (1) trapping full implicit calculations, second argument: Fokker-Planck (0) or Drift Kinetic Equation (1 or 2) [1,2]
       - ia: inverse aspect ratio [1,nr]
       - xftp_norm: normalization factor for the drift kinetic equation [1,nr]
       - XXDr: Fast electron radial diffusion (nhuth_ref*ap^2) [np,nmhu,nr]
       - XXVr: Fast electron radial pinch (nhuth_ref*ap) [np,nmhu,nr]
       - delta: Toroidal magnetic ripple depth [1,nr]
       - Nq: Number of toroidal magnetic coils * local safety factor [1,nr]
       - coll_mode: 0 -> Maxwellian relativistic collision operator, 1 -> High-velocity limit  [1,1], 2 -> linearizeded Belaiev-Budker , 3 -> Non relativistic Lorentz model
           (default = 2)
       - dtn: normalized integration time step (dtn > 0: fully implicit time differencing scheme, dtn < 0: Crank-Nicholson time differencing scheme) (1/nhu_th) [1,1]
           (default = 1000)
       - prec0: required accuracy on the residu for f and for g (%) [1,2]
           (default = [1e-10,1e-10])
       - nit: maximum number of iterations  for f and g [1,2]
           (default = [20,20])
       - invparam: structure for parameters of the matrix inversion
       - display: display mode parameter (full: 2, partial: 1, no: = 0) index of the radial position where results are displayed
           (default = [0,1])
       - XXfinit: initial guess for the 2D electron momentum distribution function [np,nmhu,nr]
           (default = Maxwellian)

   Output:
       - XXf0: 2D electron momentum distribution function normalized to the normalized density [np,nmhu,nr]
       - XXf0_tp: First neoclassical correction (f_tp) [np,nmhu,nr]
       - XXf0_g: Second neoclassical correction (g) [np,nmhu,nr]
       - XXSp0: p-component of the 2D electron momentum flux for f0 [np,nmhu,nr]
       - XXSmhu0: mhu-component of the 2D electron momentum flux for f0 [np,nmhu,nr]
       - XXSp0_g: p-component of the 2D electron momentum flux for g [np,nmhu,nr]
       - XXSmhu0_g: mhu-component of the 2D electron momentum flux for g [np,nmhu,nr]
       - XXSp0_tp: p-component of the 2D electron momentum flux for ftp [np,nmhu,nr]
       - XXSmhu0_tp: mhu-component of the 2D electron momentum flux for ftp [np,nmhu,nr]
       - Pcoll_tot: Power loss on collisions (me*ne*nhuth*vth^2) [1,nr]
       - Pcoll_tot_fsav: FS averaged power loss on collisions from f0+g+f_tp(me*ne*nhuth*vth^2) [1,nr]
       - Pabs0_fsav: FS averaged  RF absorbed power from f0(me*ne*nhuth*vth^2) [1,nr]
       - Pabs_tot_fsav: FS averaged  RF absorbed power from f0+g+f_tp(me*ne*nhuth*vth^2) [1,nr]
       - curr_0: RF + ohmic current density [1,nr]
       - curr_tot: Total current density (RF + ohmic + boostrap current) [1,nr]
       - curr_0_fsav:  FS averaged  RF + ohmic current density [1,nr]
       - curr_tot_fsav:  FS averaged  Total current density (RF + ohmic + boostrap current) [1,nr]
       - currB: Bulk current (RF + ohmic + boostrap - runway current) [1,nr]
       - RR: Runaway rate through the surface of a sphere of radius p [nr,np]
       - xepsi: normalized ohmic electric field at Bmin (pth_ref*nhuth_ref/e) [1,nr]
   - trapped_fraction: exact fraction of trapped electrons ft [1,nr]
       - Ten_out: New temperature of the bulk (mec2) [1,nr]
       - normf_tot_fsav: Flux surface averaged norm of the distribution function [1,nr]
       - curr0_fsav: Flux surface averaged RF + ohmic current density [1,nr]
       - curr_tot_fsav: Total flux surface averaged current density (RF + ohmic + boostrap current)
       - Xpn: Momentum grid normalized to the reference momentum value [np,nmhu]
       - Xmhu: Angular grid [np,nmhu]
       - mhubounce2: Square of the bounce averaged cosine angle [1,1]
       - memoryLU: memory size of the L+U matrices for f0 in MBytes[1,1]
  
   Warning: none


Fokker-Planck part originally done by Y. Peysson 
Drift kinetic calculations for arbitrary magnetic equilibrium done by J. Decker (MIT-RLE) <jodecker@mit.edu> and Y. Peysson (CEA-DRFC) <yves.peysson@cea.fr>
PETSc implementation by Frederic Teisseire (CEA-DRFC) <frederic.tesseire@cea.fr>

0001 function [XXf0,XXf0_tp,XXf0_g,XXf0_L_tp,XXf0_L_g,radial_mode,memoryLU_f,memoryLU_g,time_DKE,normf0,curr0,residu_f] = solver_dke_yp(dkepath,dkeparam,display,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,XXSavalanches,xrr_out,XXfM,XXfinit)
0002 %
0003 %
0004 %   Self-consistent 3-D relativistic solver of the drift kinetic equation in momentum space with Beliaev-Budker collision integrals
0005 %   Includes toroidal magnetic losses for the Fokker-Planck equation only for circular plasma cross-section (TS)
0006 %   Fully implicit or Crank-Nicholson time differencing schemes
0007 %   LH and EC wave quasilinear diffusion calculations + ohmic electric field
0008 %   IMPORTANT: All input quantities are defined on the distribution function grid: half grid
0009 %
0010 %   Version 4: non-uniform momentum and pitch-angle grids + fully implicit 3D calculations
0011 %
0012 %   See references:
0013 %
0014 %       - C. Karney, Comput. Phys. Rep. 4 (1986) 183.
0015 %       - B. Braams, C. Karney, Phys.fluids B1 (1989) 1355
0016 %       - M. Shoucri and I. Shkarofsky, Comp. Phys. Comm. 82 (1994) 287
0017 %       - I. Shkarofsky and M. Shoucri, Nuclear Fusion, 37 (1997) 539
0018 %       - Y. Peysson and M. Shoucri, Comp. Phys. Comm. 109 (1998) 55
0019 %       - S.D. Schultz, A. Bers and A.B. Rahm, AIP Proc. of the 12th Conf. Radio Frequency Power in Plasmas, Savannah (USA), 403 1997,327
0020 %       - M. Ju et al., Phys. Plasmas, 9 (2002) 493
0021 %
0022 %   (see test_dke_1yp.m)
0023 %
0024 %   Input:
0025 %
0026 %       - betath: Normalized reference thermal velocity (vth/c) [1,1]
0027 %       - Rp_norm: Normalized plasma major radius (ap) [1,1]
0028 %       - xTe_norm: normalized electron temperature (Te/Te_ref) [1,nr]
0029 %       - xne_norm: normalized electron density (ne/ne_ref) [1,nr]
0030 %       - Zi: ion charges [1,m]
0031 %       - mi: ion masses (u.m.a.) [1,m]
0032 %       - xni_norm: normalized ion density (ni/ne_ref) [m,nr]
0033 %       - pn: normalized electron momentum space grid (pth_ref) [1,np]
0034 %       - mhu: cos(theta) angular grid [1,nmhu]
0035 %       - XXDlh: LH normalized quasi-linear diffusion coefficient (nhuth_ref*pth_ref^2) [np,nmhu,nr]
0036 %       - XXDcy: EC normalized quasi-linear diffusion coefficient (nhuth_ref*pth_ref^2) [np,nmhu,nr]
0037 %       - param_cy: *******************normalized EC frequency [1,nr]
0038 %       - XXsinksource: normalized source/loss frequency (nhuth_ref) [np,nmhu,nr]
0039 %       - xepsi_fsav: flux surface averaged normalized ohmic electric field (pth_ref*nhuth_ref/e) [1,nr]
0040 %       - rdke: list of radial positions indexes at which FP or DKE are solved [1,nr_dke].
0041 %       - bounce: First argument: bounce averaged option, (0) no trapping, (1) trapping full implicit calculations, second argument: Fokker-Planck (0) or Drift Kinetic Equation (1 or 2) [1,2]
0042 %       - ia: inverse aspect ratio [1,nr]
0043 %       - xftp_norm: normalization factor for the drift kinetic equation [1,nr]
0044 %       - XXDr: Fast electron radial diffusion (nhuth_ref*ap^2) [np,nmhu,nr]
0045 %       - XXVr: Fast electron radial pinch (nhuth_ref*ap) [np,nmhu,nr]
0046 %       - delta: Toroidal magnetic ripple depth [1,nr]
0047 %       - Nq: Number of toroidal magnetic coils * local safety factor [1,nr]
0048 %       - coll_mode: 0 -> Maxwellian relativistic collision operator, 1 -> High-velocity limit  [1,1], 2 -> linearizeded Belaiev-Budker , 3 -> Non relativistic Lorentz model
0049 %           (default = 2)
0050 %       - dtn: normalized integration time step (dtn > 0: fully implicit time differencing scheme, dtn < 0: Crank-Nicholson time differencing scheme) (1/nhu_th) [1,1]
0051 %           (default = 1000)
0052 %       - prec0: required accuracy on the residu for f and for g (%) [1,2]
0053 %           (default = [1e-10,1e-10])
0054 %       - nit: maximum number of iterations  for f and g [1,2]
0055 %           (default = [20,20])
0056 %       - invparam: structure for parameters of the matrix inversion
0057 %       - display: display mode parameter (full: 2, partial: 1, no: = 0) index of the radial position where results are displayed
0058 %           (default = [0,1])
0059 %       - XXfinit: initial guess for the 2D electron momentum distribution function [np,nmhu,nr]
0060 %           (default = Maxwellian)
0061 %
0062 %   Output:
0063 %       - XXf0: 2D electron momentum distribution function normalized to the normalized density [np,nmhu,nr]
0064 %       - XXf0_tp: First neoclassical correction (f_tp) [np,nmhu,nr]
0065 %       - XXf0_g: Second neoclassical correction (g) [np,nmhu,nr]
0066 %       - XXSp0: p-component of the 2D electron momentum flux for f0 [np,nmhu,nr]
0067 %       - XXSmhu0: mhu-component of the 2D electron momentum flux for f0 [np,nmhu,nr]
0068 %       - XXSp0_g: p-component of the 2D electron momentum flux for g [np,nmhu,nr]
0069 %       - XXSmhu0_g: mhu-component of the 2D electron momentum flux for g [np,nmhu,nr]
0070 %       - XXSp0_tp: p-component of the 2D electron momentum flux for ftp [np,nmhu,nr]
0071 %       - XXSmhu0_tp: mhu-component of the 2D electron momentum flux for ftp [np,nmhu,nr]
0072 %       - Pcoll_tot: Power loss on collisions (me*ne*nhuth*vth^2) [1,nr]
0073 %       - Pcoll_tot_fsav: FS averaged power loss on collisions from f0+g+f_tp(me*ne*nhuth*vth^2) [1,nr]
0074 %       - Pabs0_fsav: FS averaged  RF absorbed power from f0(me*ne*nhuth*vth^2) [1,nr]
0075 %       - Pabs_tot_fsav: FS averaged  RF absorbed power from f0+g+f_tp(me*ne*nhuth*vth^2) [1,nr]
0076 %       - curr_0: RF + ohmic current density [1,nr]
0077 %       - curr_tot: Total current density (RF + ohmic + boostrap current) [1,nr]
0078 %       - curr_0_fsav:  FS averaged  RF + ohmic current density [1,nr]
0079 %       - curr_tot_fsav:  FS averaged  Total current density (RF + ohmic + boostrap current) [1,nr]
0080 %       - currB: Bulk current (RF + ohmic + boostrap - runway current) [1,nr]
0081 %       - RR: Runaway rate through the surface of a sphere of radius p [nr,np]
0082 %       - xepsi: normalized ohmic electric field at Bmin (pth_ref*nhuth_ref/e) [1,nr]
0083 %   - trapped_fraction: exact fraction of trapped electrons ft [1,nr]
0084 %       - Ten_out: New temperature of the bulk (mec2) [1,nr]
0085 %       - normf_tot_fsav: Flux surface averaged norm of the distribution function [1,nr]
0086 %       - curr0_fsav: Flux surface averaged RF + ohmic current density [1,nr]
0087 %       - curr_tot_fsav: Total flux surface averaged current density (RF + ohmic + boostrap current)
0088 %       - Xpn: Momentum grid normalized to the reference momentum value [np,nmhu]
0089 %       - Xmhu: Angular grid [np,nmhu]
0090 %       - mhubounce2: Square of the bounce averaged cosine angle [1,1]
0091 %       - memoryLU: memory size of the L+U matrices for f0 in MBytes[1,1]
0092 %
0093 %   Warning: none
0094 %
0095 %
0096 %Fokker-Planck part originally done by Y. Peysson
0097 %Drift kinetic calculations for arbitrary magnetic equilibrium done by J. Decker (MIT-RLE) <jodecker@mit.edu> and Y. Peysson (CEA-DRFC) <yves.peysson@cea.fr>
0098 %PETSc implementation by Frederic Teisseire (CEA-DRFC) <frederic.tesseire@cea.fr>
0099 %
0100 if nargin < 9,
0101     error('Wrong number of input arguments in solver_dke_yp')
0102 end
0103 %
0104 time00 = clock;
0105 timeCPU0 = cputime;
0106 memoryLU_f = 0;
0107 memoryLU_g = 0;
0108 %
0109 boundary_mode_f = dkeparam.boundary_mode_f;
0110 boundary_mode_g = dkeparam.boundary_mode_g;
0111 dke_mode = dkeparam.dke_mode;
0112 bounce_mode = dkeparam.bounce_mode;
0113 dtn = dkeparam.dtn;
0114 temppath = dkepath.temppath;
0115 prec0_f = dkeparam.prec0_f;
0116 prec0_g = dkeparam.prec0_g;
0117 nit_f = dkeparam.nit_f;
0118 nit_g = dkeparam.nit_g;
0119 %
0120 display_mode = display.display_mode;
0121 %
0122 nr = gridDKE.nr;
0123 npn = gridDKE.npn;
0124 nmhu = gridDKE.nmhu;
0125 nr_dke = gridDKE.nr_dke;
0126 mhu = gridDKE.mhu;
0127 Xdmhumm = gridDKE.Xdmhumm;
0128 Xdmhum = gridDKE.Xdmhum;
0129 Xdmhu = gridDKE.Xdmhu;
0130 Xdmhup = gridDKE.Xdmhup;
0131 Xdmhupp = gridDKE.Xdmhupp;
0132 Xdpnmm = gridDKE.Xdpnmm;
0133 Xdpnm = gridDKE.Xdpnm;
0134 Xdpn = gridDKE.Xdpn;
0135 Xdpnp = gridDKE.Xdpnp;
0136 Xdpnpp = gridDKE.Xdpnpp;
0137 Xmhumm = gridDKE.Xmhumm;
0138 Xmhum = gridDKE.Xmhum;
0139 Xmhu = gridDKE.Xmhu;
0140 Xmhup = gridDKE.Xmhup;
0141 Xmhupp = gridDKE.Xmhupp;
0142 X1mmhu2m = gridDKE.X1mmhu2m;
0143 X1mmhu2 = gridDKE.X1mmhu2;
0144 X1mmhu2p = gridDKE.X1mmhu2p;
0145 Xpnmm = gridDKE.Xpnmm;
0146 Xpnm = gridDKE.Xpnm;
0147 Xpn = gridDKE.Xpn;
0148 Xpnp = gridDKE.Xpnp;
0149 Xpnpp = gridDKE.Xpnpp;
0150 Xpn2m = gridDKE.Xpn2m;
0151 Xpn2 = gridDKE.Xpn2;
0152 XXpn2 = gridDKE.XXpn2;
0153 Xpn2p = gridDKE.Xpn2p;
0154 xpsin = gridDKE.xpsin;
0155 rdke = gridDKE.rdke;
0156 xrho = gridDKE.xrho;
0157 %
0158 pn = gridDKE.pn;
0159 pn2 = gridDKE.pn2;
0160 pn3 = gridDKE.pn3;
0161 xpn = gridDKE.xpn;
0162 xpn2 = gridDKE.xpn2;
0163 xpn3 = gridDKE.xpn3;dpn = gridDKE.dpn;
0164 mhu = gridDKE.mhu;
0165 dmhu = gridDKE.dmhu;
0166 pn2 = gridDKE.pn2;
0167 masku = gridDKE.masku;
0168 maskl = gridDKE.maskl;
0169 ir_display = gridDKE.ir_display_out;
0170 %
0171 z = Zmomcoef.z;
0172 z2 = Zmomcoef.z2;
0173 xz = Zmomcoef.xz;
0174 xz2 = Zmomcoef.xz2;
0175 gamma = Zmomcoef.gamma;
0176 xgamma = Zmomcoef.xgamma;
0177 Xgamma = Zmomcoef.Xgamma;
0178 sigma = Zmomcoef.sigma;
0179 xsigma = Zmomcoef.xsigma;
0180 J1 = Zmomcoef.J1;
0181 J2 = Zmomcoef.J2;
0182 J3 = Zmomcoef.J3;
0183 J4 = Zmomcoef.J4;
0184 xJ1 = Zmomcoef.xJ1;
0185 xJ2 = Zmomcoef.xJ2;
0186 xJ3 = Zmomcoef.xJ3;
0187 xJ4 = Zmomcoef.xJ4;
0188 %
0189 XXheaviside = Zbouncecoef.XXheaviside;
0190 XXRlambda_b_p1p1 = Zbouncecoef.XXRlambda_b_p1p1;
0191 XXlambda_b_p2 = Zbouncecoef.XXlambda_b_p2;
0192 XXlambda = Zbouncecoef.XXlambda;
0193 Zmask_f0_t = Zbouncecoef.Zmask_f0_t;
0194 Zmask_f0_Tm = Zbouncecoef.Zmask_f0_Tm;
0195 Zmask_f0_Tp = Zbouncecoef.Zmask_f0_Tp;
0196 ns_f = Zbouncecoef.ns_f; 
0197 Zmask_g = Zbouncecoef.Zmask_g;
0198 ns_g = Zbouncecoef.ns_g;
0199 xqtilde = Zbouncecoef.xqtilde;
0200 xqhat =  Zbouncecoef.xqhat;
0201 xq = Zbouncecoef.xq;
0202 xqbar =  Zbouncecoef.xqbar;
0203 %
0204 MMXR_t = [];
0205 %
0206 xB0 = equilDKE.xB0;               
0207 %
0208 betath_ref = mksa.betath_ref;
0209 xftp_norm = mksa.xftp_norm_ref;
0210 %
0211 [xTe_norm,xne_norm,xzni_norm,xzTi_norm,xnloss_norm,xbetath,xlnc_e,xnhu,xrnhuth] = normcoef_dke_yp(mksa,equilDKE,gridDKE);   
0212 %
0213 if dtn > 0,
0214     time_scheme = 1;%Implicit differencing scheme
0215 elseif dkeparam.dtn < 0,
0216     dtn  = abs(dtn);
0217     time_scheme = 2;%Crank-Nicholson differencing scheme
0218 end 
0219 %
0220 if isempty(MMXR_f_t) == 0, radial_mode = 1;else radial_mode = 0;end
0221 if ~isfield(dkeparam,'opsplit'),dkeparam.opsplit = 0;end
0222 %
0223 spparms('default');%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions only
0224 spparms('supernd',4);%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions only
0225 spparms('rreduce',4);%Set parameters for sparse matrix ordering algorithms and direct linear equation solver MatLab build-in functions only
0226 %
0227 if isempty(which('initmumps')),
0228     dkeparam.invproc = 1;%MatLab CGS iterative matrix inversion method enforced
0229     disp('**************************************************************************************************************');
0230     disp('WARNING: MUMPS solver called but not installed for LUKE. MatLab CGS iterative matrix inversion method enforced');
0231     disp('**************************************************************************************************************');
0232 end
0233 %
0234 %************************************************** Simulation parameter display ********************************************************
0235 %
0236 if display_mode >=1,
0237     disp('-------------------------------------------------------------------------------------------------------------------');
0238     disp('-------------------------------------------------------------------------------------------------------------------');
0239     disp('Kinetic calculation parameters: ');
0240 %Matrix facorization technique
0241     if (dkeparam.invproc > 0),        disp('   - MatLab incomplete LU matrix factorization');
0242     elseif (dkeparam.invproc ==- 2),  disp('   - MUMPSMEX incomplete LU matrix factorization');        
0243     elseif (dkeparam.invproc < -1 & dkeparam.invproc ~=- 2),  disp('   - Direct matrix inversion by external solvers');end
0244 %Operator splitting
0245     if (dkeparam.opsplit == 1),       disp('   - Operator splitting between momentum and radial dynamics.');
0246     elseif (dkeparam.opsplit == 0),       disp('   - No operator splitting between momentum and radial dynamics.');end
0247 %Numerical matrix inversion procedure
0248     if (dkeparam.invproc == 1),  disp('   - MatLab CGS iterative matrix inversion method');
0249     elseif (dkeparam.invproc == 2),  disp('   - MatLab BICG iterative matrix inversion method');
0250     elseif (dkeparam.invproc == 3),  disp('   - MatLab QMR iterative matrix inversion method');
0251     elseif (dkeparam.invproc == -1),  disp('   - MUMPS matrix inversion - FORTRAN 95 parallel processing');
0252     elseif (dkeparam.invproc == -2),  disp('   - MUMPSMEX matrix inversion - FORTRAN 95 sequential processing.');
0253     elseif (dkeparam.invproc == -3),  disp('   - SUPERLUMEX matrix inversion - C sequential processing.');
0254     elseif (dkeparam.invproc == -4),  disp('   - PETSc matrix inversion. Several solvers may be used in this package.');end
0255 %Boundary and normalization conditions
0256     if real(dkeparam.norm_mode_f) == 0,disp('   - No normalization of f0 at each iteration and radius: fully natural conservative mode');
0257     elseif real(dkeparam.norm_mode_f) == 1,    disp('   - Normalization of f0 with respect to reference density at each iteration and radius');
0258     elseif real(dkeparam.norm_mode_f) == 2,disp('   - No normalization of f0 with respect to initial density at each iteration and radius: fully natural conservative mode');end
0259     if imag(dkeparam.norm_mode_f) == 1,    disp('   - Normalization of f0 accounting for boundary fluxes at each iteration and radius');end
0260     if (dkeparam.boundary_mode_f >= 1),    disp(['   - Maxwellian distribution for FPE enforced on ',int2str(dkeparam.boundary_mode_f),' momentum grid point around p = 0']);
0261     elseif (dkeparam.boundary_mode_f == 0),disp('   - No Maxwellian distribution for FPE enforced around p = 0');end
0262     if (dkeparam.boundary_mode_g >= 1) & (dke_mode == 1),    disp(['   - Lorentz model distribution for DKE enforced on ',int2str(dkeparam.boundary_mode_g),' momentum grid point around p = 0']);
0263     elseif (dkeparam.boundary_mode_g == 0) & (dke_mode == 1),disp('   - No Lorentz model distribution for DKE enforced around p = 0');end
0264 %Time scheme
0265     if (time_scheme == 1),    disp(['   - Implicit time differencing scheme']);
0266     elseif (time_scheme == 2),disp(['   - Crank-Nicholson time differencing scheme']);end
0267     disp(['   - Normalized time step: ',int2str(dkeparam.dtn)]);
0268 %Momentum dynamics
0269     if (dkeparam.coll_mode == 0),     disp(['   - Collisions: Relativistic Maxwellian background']);
0270     elseif (dkeparam.coll_mode == 1), disp(['   - Collisions: High velocity limit']);    
0271     elseif (dkeparam.coll_mode == 2), disp(['   - Collisions: Linearized relativistic Belaiev-Budker']);
0272     elseif (dkeparam.coll_mode == 3), disp(['   - Collisions: Non-relativistic Lorentz model']);
0273     elseif (dkeparam.coll_mode == 4), disp(['   - Collisions: Non-relativistic Maxwellian background']);end
0274     if (bounce_mode == 0),      disp(['   - No bounce averaging']);
0275     elseif (bounce_mode == 1),  disp(['   - Bounce averaging + implicit symmetrization in the bounce domain (15 diagonals operator) for f0']);end
0276     if (dke_mode == 0),       disp(['   - Resolution of the electron Fokker-Planck Equation']);
0277     elseif (dke_mode == 1),   disp(['   - Resolution of the electron Drift Kinetic Equation']);end
0278 %Spatial dynamics
0279     if (radial_mode == 1),      disp(['   - Fast electron radial diffusion or pinch']);
0280     elseif (radial_mode == 0),  disp(['   - No fast electron radial diffusion or pinch']);end
0281 %
0282     disp(['   - Momentum grid: ',int2str(npn)]);
0283     disp(['   - Angular grid: ',int2str(nmhu)]);
0284     disp(['   - Total radial grid: ',int2str(nr)]);
0285     disp(['   - Effective radial grid: ',int2str(nr_dke)]);
0286     disp('-------------------------------------------------------------------------------------------------------------------');
0287 end
0288 %
0289 %normf0 = zeros(nr,nit_f);
0290 normf1 = zeros(nr,nit_f);
0291 %curr0 = zeros(nr,nit_f);
0292 curr1 = zeros(nr,nit_f);
0293 prec_curr = zeros(nr,nit_f);
0294 prec_normf = zeros(nr,nit_f);
0295 %residu_f = NaN*ones(nr,nit_f);
0296 prec_residu_f = NaN*ones(nr,nit_f);
0297 residu_g = NaN*ones(nr_dke,nit_g);
0298 prec_residu_g = NaN*ones(nr_dke,nit_g);
0299 %
0300 XXI = zeros(npn,nmhu,nr);
0301 XXI_g = zeros(npn,nmhu,nr);
0302 XXf0 = zeros(npn,nmhu,nr);
0303 XXf0_tp = zeros(npn,nmhu,nr);
0304 XXf0_g = zeros(npn,nmhu,nr);
0305 XXf0_L_tp = zeros(npn,nmhu,nr);
0306 XXf0_L_g = zeros(npn,nmhu,nr);  
0307 XXf0_g_th = zeros(npn,nmhu,nr);
0308 %
0309 if ~isfield(dkeparam.PETScparam,'cloop'),%Backward compatibility
0310     warning('dkeparam.PETScparam.cloop field does not exists. It is created and the calculation is done with Matlab Built-in function. Please fixe the missing field in the dkeparam structure');    
0311     dkeparam.PETScparam.cloop = 0;
0312 end
0313 %
0314 time0 = clock;
0315 %
0316 if dkeparam.opsplit == 0,
0317     MMXdiag = diag(MMXP_f_t/time_scheme + MMXR_f_t/time_scheme + MMXT_f_t);
0318 %
0319     if isempty(find(MMXdiag==0)) == 1, 
0320         MMXPDI_f_t = inv(diag(MMXdiag));
0321         MMXPC_f_t = MMXPDI_f_t*(MMXP_f_t/time_scheme + MMXR_f_t/time_scheme + MMXT_f_t);%Main diagonal matrix preconditionning
0322 %
0323         if dkeparam.invproc >=1,%MatLab build-in matrix incomplete factorization only
0324             [MMXL_f_t,MMXU_f_t] = ilu(MMXPC_f_t,struct('type','crout','milu','row','droptol',dkeparam.ludroptol));
0325         elseif dkeparam.invproc == -2,%MUMPSMEX case only (single processor). See MUMPS documentation
0326             job_MUMPSMEX = 1;MUMPSMEX_dke_yp;
0327             job_MUMPSMEX = 2;MUMPSMEX_dke_yp;
0328         end
0329 %
0330         time_f = etime(clock,time0);
0331         if dkeparam.invproc >=1,
0332             memoryLU_f = (getfield(whos('MMXL_f_t'),'bytes') + getfield(whos('MMXU_f_t'),'bytes'))/1000000;%MBytes
0333         elseif dkeparam.invproc == -2,
0334             memoryLU_f = (id_MUMPSMEX.INFOG(9)*8 + id_MUMPSMEX.INFOG(10)*4)/1000000;%MBytes, see MUMPS documentation (8 bytes for double precision float, 4 byte for integer)
0335         end
0336     else
0337         error('Fokker-Planck matrix has a null element on the main diagonal');%Matrix must be main diagonal dominant
0338     end
0339 %
0340     if dke_mode == 1,
0341         time0 = clock;  
0342         MMXPDI_g_t = inv(diag(diag(MMXP_g_t)));
0343         MMXPC_g_t = MMXPDI_g_t*MMXP_g_t;%Main diagonal matrix preconditionning
0344     %
0345         if dkeparam.invproc >=1,%MatLab build-in matrix incomplete factorization only
0346            [MMXL_g_t,MMXU_g_t] = ilu(MMXPC_g_t,struct('type','crout','milu','row','droptol',dkeparam.ludroptol)); 
0347         end
0348     %
0349         time_g = etime(clock,time0);
0350         if dkeparam.invproc >=1,
0351            memoryLU_g = (getfield(whos('MMXL_g_t'),'bytes') + getfield(whos('MMXU_g_t'),'bytes'))/1000000;
0352         elseif dkeparam.invproc == -2,
0353            memoryLU_g = (id_MUMPSMEX.INFOG(9)*8 + id_MUMPSMEX.INFOG(10)*4)/1000000;%MBytes, see MUMPS documentation (8 bytes for double precision float, 4 byte for integer)
0354         end
0355     end
0356     %
0357     if display_mode >= 1,   
0358         if dkeparam.invproc >=1,
0359             info_dke_yp(2,['Size of the MatLab LU matrices for f: ',num2str(memoryLU_f),' MBytes']);
0360             info_dke_yp(2,['Full elapsed time of MatLab LU matrix factorization for f: ',num2str(time_f),' (s)']);
0361         elseif dkeparam.invproc == -2,
0362             info_dke_yp(2,'LU matrix factorization within MUMPSMEX');
0363         elseif ~dkeparam.invproc == -2,
0364             info_dke_yp(2,'Direct matrix inversion, no LU matrix factorization');
0365         end
0366         if dke_mode == 1,
0367             if dkeparam.invproc >=1,
0368                 info_dke_yp(2,['Size of the MatLab LU matrices for g: ',num2str(memoryLU_g),' MBytes']);
0369                 info_dke_yp(2,['Full elapsed time of MatLab LU matrix factorization for g: ',num2str(time_g),' (s)']);
0370             elseif dkeparam.invproc == -2,
0371                 info_dke_yp(2,'LU matrix factorization within MUMPSMEX');
0372             elseif ~dkeparam.invproc == -2,
0373                 info_dke_yp(2,'Direct matrix inversion, no LU matrix factorization');
0374             end
0375         end
0376     end
0377 elseif dkeparam.opsplit == 1,
0378     MMXPdiag = diag(MMXP_f_t/time_scheme + 2*MMXT_f_t);
0379     MMXRdiag = diag(MMXR_f_t/time_scheme + 2*MMXT_f_t);
0380     %
0381     if isempty(find(MMXPdiag==0)) == 1 | isempty(find(MMXRdiag==0)) == 1, 
0382         MMXPDI_P_f_t = inv(diag(MMXPdiag));
0383         MMXPDI_R_f_t = inv(diag(MMXRdiag));
0384         %
0385         MMXPC_P_f_t = MMXPDI_P_f_t*(MMXP_f_t/time_scheme + 2*MMXT_f_t);%Main diagonal matrix preconditionning
0386         MMXPC_R_f_t = MMXPDI_R_f_t*(MMXR_f_t/time_scheme + 2*MMXT_f_t);%Main diagonal matrix preconditionning
0387 %
0388         if dkeparam.invproc >=1,%MatLab build-in matrix incomplete factorization only
0389             [MMXL_P_f_t,MMXU_P_f_t] = ilu(MMXPC_P_f_t,struct('type','crout','milu','row','droptol',dkeparam.ludroptol));
0390             [MMXL_R_f_t,MMXU_R_f_t] = ilu(MMXPC_R_f_t,struct('type','crout','milu','row','droptol',dkeparam.ludroptol));
0391         elseif dkeparam.invproc == -2,%MUMPSMEX case only (single processor). See MUMPS documentation
0392             job_MUMPSMEX = 1;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;%P matrix
0393             job_MUMPSMEX = 2;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;%P matrix
0394             %
0395             if dkeparam.invproc >=1,
0396                 memoryLU_P_f = (getfield(whos('MMXL_P_f_t'),'bytes') + getfield(whos('MMXU_P_f_t'),'bytes'))/1000000;%MBytes
0397             elseif dkeparam.invproc == -2,
0398                 memoryLU_P_f = (id_MUMPSMEX_P.INFOG(9)*8 + id_MUMPSMEX_P.INFOG(10)*4)/1000000;%MBytes, see MUMPS documentation (8 bytes for double precision float, 4 byte for integer)
0399             end
0400             %
0401             job_MUMPSMEX = 1;mode_MUMPSMEX=0;MUMPSMEX_dke_yp;%R matrix
0402             job_MUMPSMEX = 2;mode_MUMPSMEX=0;MUMPSMEX_dke_yp;%R matrix
0403             %
0404             if dkeparam.invproc >=1,
0405                 memoryLU_R_f = (getfield(whos('MMXL_R_f_t'),'bytes') + getfield(whos('MMXU_R_f_t'),'bytes'))/1000000;%MBytes
0406             elseif dkeparam.invproc == -2,
0407                 memoryLU_R_f = (id_MUMPSMEX_R.INFOG(9)*8 + id_MUMPSMEX_R.INFOG(10)*4)/1000000;%MBytes, see MUMPS documentation (8 bytes for double precision float, 4 byte for integer)
0408             end
0409         end
0410 %
0411         time_f = etime(clock,time0);
0412     else
0413         error('One of the Fokker-Planck two matrices has a null element on the main diagonal');%Matrix must be main diagonal dominant
0414     end
0415 elseif dkeparam.opsplit == 2,    
0416     MMXPdiag = diag(MMXP_f_t/time_scheme + MMXT_f_t);
0417     %
0418     if isempty(find(MMXPdiag==0)) == 1, 
0419         MMXPDI_P_f_t = inv(diag(MMXPdiag));
0420         %
0421         MMXPC_P_f_t = MMXPDI_P_f_t*(MMXP_f_t/time_scheme + MMXT_f_t);%Main diagonal matrix preconditionning
0422 %
0423         if dkeparam.invproc >=1,%MatLab build-in matrix incomplete factorization only
0424             [MMXL_P_f_t,MMXU_P_f_t] = ilu(MMXPC_P_f_t,struct('type','crout','milu','row','droptol',dkeparam.ludroptol));
0425         elseif dkeparam.invproc == -2,%MUMPSMEX case only (single processor). See MUMPS documentation
0426             job_MUMPSMEX = 1;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;%P matrix
0427             job_MUMPSMEX = 2;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;%P matrix
0428             %
0429             if dkeparam.invproc >=1,
0430                 memoryLU_P_f = (getfield(whos('MMXL_P_f_t'),'bytes') + getfield(whos('MMXU_P_f_t'),'bytes'))/1000000;%MBytes
0431             elseif dkeparam.invproc == -2,
0432                 memoryLU_P_f = (id_MUMPSMEX_P.INFOG(9)*8 + id_MUMPSMEX_P.INFOG(10)*4)/1000000;%MBytes, see MUMPS documentation (8 bytes for double precision float, 4 byte for integer)
0433             end
0434             %
0435         end
0436 %
0437         time_f = etime(clock,time0);
0438     else
0439         error('One of the Fokker-Planck two matrices has a null element on the main diagonal');%Matrix must be main diagonal dominant
0440     end
0441 end
0442 %
0443 if display_mode == 2,
0444     figure('Name','Momentum space matrix'),spy(MMXP_f_t);
0445     zoom on;
0446     axis('square');
0447     title('Momentum space matrix');
0448 %
0449     if radial_mode == 1,
0450         figure('Name','Configuration space matrix'),spy(MMXR_f_t);
0451         zoom on;
0452         axis('square');
0453         title('Configuration space matrix');    
0454     end
0455 %
0456     figure('Name','Total preconditionned matrix'),spy(MMXPC_f_t);
0457     zoom on;
0458     axis('square');
0459     title('Total preconditionned matrix');
0460 end
0461 %
0462 %************************************************** Main Fokker-Planck loop ********************************************************
0463 %
0464 if display_mode >= 1,
0465     disp('-------------------------------------------------------------------------------------------------------------------');
0466 end
0467 %
0468 time0 = clock;
0469 XXf0 = XXfinit;%Initial guess (normalized to xne_norm)
0470 %
0471 % The avalanche operator is assumed to be conservative.
0472 % The corresponding source term is compensated for by removing an equivalent fraction of the Maxwellian
0473 %
0474 xnormavalanches = xrr_out + 2*pi*(xqtilde./xqhat).*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda.*XXpn2.*XXSavalanches,1),1).';
0475 %
0476 % For differential operators, the TP symmetrization and the multiplication
0477 % by pn^2 is done in matrixcalc. For the sink/source terms, it is done here for
0478 % the sake of simplicity.
0479 %
0480 XXsinksource_pn2 = XXpn2.*(XXsinksource + XXSavalanches);
0481 %
0482 if bounce_mode > 0,
0483     for ir = 1:nr,
0484         XXsinksource_pn2(:,Zmask_f0_Tp{ir},ir) = (XXsinksource_pn2(:,Zmask_f0_Tp{ir},ir) + fliplr(XXsinksource_pn2(:,Zmask_f0_Tm{ir},ir)))/2;
0485         XXsinksource_pn2(:,Zmask_f0_Tm{ir},ir) = fliplr(XXsinksource_pn2(:,Zmask_f0_Tp{ir},ir));% to avoid confusion if writeCloopinput_dke_ft is used.
0486     end
0487 end
0488 %
0489 % Initial distribution
0490 %
0491 MMXfinit_t = [];
0492 for ir = 1:nr,
0493     MMXfinit_t = [MMXfinit_t;reshape(XXfinit(:,Zmask_f0_t{ir},ir)',(nmhu - ns_f(ir))*npn,1)];
0494 end
0495 %
0496 if dkeparam.PETScparam.cloop & dkeparam.invproc == -4,
0497     %
0498 %    writeCloopinput_dke_ft(nit_f,boundary_mode_f,bounce_mode,dkeparam.coll_mode,real(dkeparam.norm_mode_f),prec0_f,display_mode,betath_ref,xTe_norm,xne_norm,dtn,mhu,nmhu,npn,ns_f,nr,sm_f,...
0499      %   masku,maskl,dpn,dmhu,xpn,xpn2,xpn3,xsigma,xgamma,xz,xz2,pn,pn2,pn3,...
0500      %   sigma,gamma,z,z2,xJ1,xJ2,xJ3,xJ4,J1,J2,J3,J4,...
0501     %    Xmhu,XXsinksource_pn2,Xmask_r_edge,XXfM,XXfinit,XXRlambda_b_p1p1,rdke,Xpn2,Zpn2_t,...
0502     %    Zmask_f0_t,...
0503     %    MMXT_f_t,MMXR_f_t,MMXP_f_t,MMXPC_f_t,MMXPDI_f_t,...
0504     %    XXI,XXf0,time_scheme,dkepath);
0505     %
0506 %    [XXf0] = PETSc_Cloop_dke_ft(dkeparam.PETScparam,display_mode,dkepath);
0507 else
0508     %
0509     for ir = 1:nr,
0510         %
0511         normf0(ir,1) = 2*pi*(xqtilde(ir)./xqhat(ir))*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,ir).*XXf0(:,:,ir).*Xpn2,1),1);
0512         curr0(ir,1) = 2*pi*(xq(ir)./xqbar(ir))*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXf0(:,:,ir).*Xpn2.*Xpn.*Xmhu./Xgamma,1),1);
0513         %
0514     end
0515     %
0516     if display_mode == 2,
0517         Fig_FPEconvstatus = figure('Name','Fokker-Planck equation convergence status');
0518     end
0519     %
0520     for it_f = 1:nit_f,     
0521         %
0522         for ir = 1:nr,
0523             if dkeparam.coll_mode == 2,%Belaiev-Budker collision operator
0524                 f1 = 1.5*integral_dke_jd(dmhu,XXf0(:,:,ir).*(ones(npn,1)*mhu),2)';               
0525                 [I1,I2] = firstordercollop_dke_yp(f1,betath_ref,...
0526                         masku,maskl,...
0527                         dpn,xpn,xpn2,xpn3,xsigma,xgamma,xz,xz2,pn,pn2,pn3,sigma,gamma,z,z2,...
0528                         xJ1,xJ2,xJ3,xJ4,J1,J2,J3,J4,...
0529                         xTe_norm(ir));
0530                 XI = 4*pi*((XXfM(:,1,ir)'.*(f1 + I1./pn2 + I2.*pn)./gamma)'*mhu).*XXRlambda_b_p1p1(:,:,ir);   
0531                 XI(1,:) = 0;%Boundary conditions at p = 0
0532                 XXI(:,:,ir) = XI;
0533             else
0534                 XXI(:,:,ir) = zeros(npn,nmhu);
0535             end
0536         end
0537         %
0538         MMXf0_t = [];
0539         MMXq_t = [];
0540         MMXIpn2_t = [];
0541         MMXsource_t = [];
0542         %
0543         for ir = 1:nr,
0544             XI_t = XXI(:,Zmask_f0_t{ir},ir);
0545             %
0546             if boundary_mode_f > 0,
0547                 XI_t(1:boundary_mode_f,:) = 0;
0548             end
0549             %
0550             MMXf0_t = [MMXf0_t;reshape(XXf0(:,Zmask_f0_t{ir},ir).',(nmhu - ns_f(ir))*npn,1)];
0551             MMXsource_t = [MMXsource_t;reshape(XXsinksource_pn2(:,Zmask_f0_t{ir},ir).',(nmhu - ns_f(ir))*npn,1)];
0552             %
0553             if ir == nr & ~isempty(Xmask_r_edge),%Edge radial boundary condition: XI_t = 0 if radial transport operates
0554                 %TODO: justifiying edge conditions for the Legendre integral
0555                 XI_t = XI_t*0;
0556             end
0557             %
0558             MMXIpn2_t = [MMXIpn2_t;reshape((XI_t.*Zpn2_t{ir})',(nmhu - ns_f(ir))*npn,1)];%Zpn2_t is here as part of the partial Jacobian used to avoid singularities at p = 0 and psi = 0, and for pitch-angle derivative df/dmhu = 0 enforced by symmetry at mhu == pi/2
0559         end 
0560         %
0561         if dkeparam.opsplit == 0,
0562             MMXq_t = MMXsource_t + MMXIpn2_t + (MMXT_f_t - (time_scheme - 1)*(MMXP_f_t + MMXR_f_t)/time_scheme)*MMXfinit_t;%Right handside term of the Fokker-Planck equation
0563             MMXR_t = sparse(MMXPDI_f_t*MMXq_t - MMXPC_f_t*MMXfinit_t);%Residu calculation for matrix inversion, force sparse type for the RHS, optimized for Matlab
0564             %
0565             %          save testLU3D.mat MMXPC_f_t MMXR_t dkeparam  MMXL_f_t MMXU_f_t dkepath display_mode
0566             %
0567             if dkeparam.invproc == 0,%Matlab direct inversion
0568                 MMXS_t = (MMXR_t\MMXPC_f_t).';
0569                 flaginv_f = 0;
0570             elseif dkeparam.invproc == 1,%Matlab cgs iterative inversion method
0571                 [MMXS_t,flaginv_f] = cgs(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t);
0572             elseif dkeparam.invproc == 2,%Matlab bicg iterative inversion method
0573                 [MMXS_t,flaginv_f] = bicg(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t); 
0574             elseif dkeparam.invproc == 3,%Matlab qmr iterative inversion method
0575                 [MMXS_t,flaginv_f] = qmr(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t);
0576             elseif dkeparam.invproc == -1,%MUMPS FORTAN 95 parallel or sequential processing package (direct inversion)
0577                 [MMXS_t,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_f_t,full(MMXR_t),display_mode,dkepath);
0578             elseif dkeparam.invproc == -2,%MUMPS FORTAN 95 parallel processing package (inversion in MEX file, with LU matrix preconditionning)
0579                 job_MUMPSMEX = 3;MUMPSMEX_dke_yp;
0580                 MMXS_t = id_MUMPSMEX.SOL; 
0581                 flaginv_f = flag_MUMPSMEX; 
0582             elseif dkeparam.invproc == -3,%SUPERLU 3.0 sequential MEX file (direct inversion)
0583                 [MMXS_t] = lusolve(MMXPC_f_t,MMXR_t);flaginv_f = 0;
0584             elseif dkeparam.invproc == -4,%Solving with PETSc
0585                 [MMXS_t,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),display_mode,dkepath);
0586             end
0587         elseif dkeparam.opsplit == 1,        
0588             MMXq_t = MMXsource_t + MMXIpn2_t + (2*MMXT_f_t - MMXR_f_t/time_scheme - (time_scheme - 1)*MMXP_f_t/time_scheme)*MMXfinit_t;%Right handside term of the Fokker-Planck equation
0589             MMXR_t = sparse(MMXPDI_P_f_t*MMXq_t - MMXPC_P_f_t*MMXfinit_t);%Residu calculation for matrix inversion, force sparse type for the RHS, optimized for Matlab
0590             %
0591             if dkeparam.invproc == 0,%Matlab direct inversion
0592                 MMXS_t = (MMXR_t\MMXPC_P_f_t).';
0593             elseif dkeparam.invproc == 1,%Matlab cgs iterative inversion method
0594                 [MMXS_t,flaginv_f] = cgs(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t);
0595             elseif dkeparam.invproc == 2,%Matlab bicg iterative inversion method
0596                 [MMXS_t,flaginv_f] = bicg(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t); 
0597             elseif dkeparam.invproc == 3,%Matlab qmr iterative inversion method
0598                 [MMXS_t,flaginv_f] = qmr(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t);
0599             elseif dkeparam.invproc == -1,%MUMPS FORTAN 95 parallel or sequential processing package (direct inversion)
0600                 [MMXS_t,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_P_f_t,full(MMXR_t),display_mode,dkepath);
0601             elseif dkeparam.invproc == -2,%MUMPS FORTAN 95 parallel processing package (inversion in MEX file, with LU matrix preconditionning)
0602                 job_MUMPSMEX = 3;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;
0603                 MMXS_t = id_MUMPSMEX_P.SOL; 
0604                 flaginv_f = flag_MUMPSMEX; 
0605             elseif dkeparam.invproc == -3,%SUPERLU 3.0 sequential MEX file (direct inversion)
0606                 [MMXS_t] = lusolve(MMXPC_P_f_t,MMXR_t);flaginv_f = 0;
0607             elseif dkeparam.invproc == -4,%Solving with PETSc
0608                 [MMXS_t,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_P_f_t,full(MMXR_t),display_mode,dkepath);
0609             end
0610             %
0611             MMXf0_t = MMXfinit_t + MMXS_t;
0612             %
0613             MMXq_t = MMXsource_t + MMXIpn2_t + (2*MMXT_f_t - MMXP_f_t/time_scheme - (time_scheme - 1)*MMXR_f_t/time_scheme)*MMXf0_t;%Right handside term of the Fokker-Planck equation
0614             MMXR_t = sparse(MMXPDI_R_f_t*MMXq_t - MMXPC_R_f_t*MMXfinit_t);%Residu calculation for matrix inversion, force sparse type for the RHS, optimized for Matlab
0615             %
0616             if dkeparam.invproc == 0,%Matlab direct inversion
0617                 MMXS_t = (MMXR_t\MMXPC_R_f_t).';
0618                 flaginv_f = 0;
0619             elseif dkeparam.invproc == 1,%Matlab cgs iterative inversion method
0620                 [MMXS_t,flaginv_f] = cgs(MMXPC_R_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_R_f_t,MMXU_R_f_t);
0621             elseif dkeparam.invproc == 2,%Matlab bicg iterative inversion method
0622                 [MMXS_t,flaginv_f] = bicg(MMXPC_R_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_R_f_t,MMXU_R_f_t); 
0623             elseif dkeparam.invproc == 3,%Matlab qmr iterative inversion method
0624                 [MMXS_t,flaginv_f] = qmr(MMXPC_R_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_R_f_t,MMXU_R_f_t);
0625             elseif dkeparam.invproc == -1,%MUMPS FORTAN 95 parallel or sequential processing package (direct inversion)
0626                 [MMXS_t,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_R_f_t,full(MMXR_t),display_mode,dkepath);
0627             elseif dkeparam.invproc == -2,%MUMPS FORTAN 95 parallel processing package (inversion in MEX file, with LU matrix preconditionning)
0628                 job_MUMPSMEX = 3;mode_MUMPSMEX=0;MUMPSMEX_dke_yp;
0629                 MMXS_t = id_MUMPSMEX_R.SOL; 
0630                 flaginv_f = flag_MUMPSMEX; 
0631             elseif dkeparam.invproc == -3,%SUPERLU 3.0 sequential MEX file (direct inversion)
0632                 [MMXS_t] = lusolve(MMXPC_R_f_t,MMXR_t);flaginv_f = 0;
0633             elseif dkeparam.invproc == -4,%Solving with PETSc
0634                 [MMXS_t,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_R_f_t,full(MMXR_t),display_mode,dkepath);
0635             end   
0636         elseif dkeparam.opsplit == 2, 
0637             MMXq_t = MMXsource_t + MMXIpn2_t + (MMXT_f_t - MMXR_f_t/time_scheme - (time_scheme - 1)*MMXP_f_t/time_scheme)*MMXfinit_t;%Right handside term of the Fokker-Planck equation
0638             MMXR_t = sparse(MMXPDI_P_f_t*MMXq_t - MMXPC_P_f_t*MMXfinit_t);%Residu calculation for matrix inversion, force sparse type for the RHS, optimized for Matlab
0639             %
0640             if dkeparam.invproc == 0,%Matlab direct inversion
0641                 MMXS_t = (MMXR_t\MMXPC_P_f_t).';
0642                 flaginv_f = 0;
0643             elseif dkeparam.invproc == 1,%Matlab cgs iterative inversion method
0644                 [MMXS_t,flaginv_f] = cgs(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t);
0645             elseif dkeparam.invproc == 2,%Matlab bicg iterative inversion method
0646                 [MMXS_t,flaginv_f] = bicg(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t); 
0647             elseif dkeparam.invproc == 3,%Matlab qmr iterative inversion method
0648                 [MMXS_t,flaginv_f] = qmr(MMXPC_P_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_P_f_t,MMXU_P_f_t);
0649             elseif dkeparam.invproc == -1,%MUMPS FORTAN 95 parallel or sequential processing package (direct inversion)
0650                 [MMXS_t,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_P_f_t,full(MMXR_t),display_mode,dkepath);
0651             elseif dkeparam.invproc == -2,%MUMPS FORTAN 95 parallel processing package (inversion in MEX file, with LU matrix preconditionning)
0652                 job_MUMPSMEX = 3;mode_MUMPSMEX=1;MUMPSMEX_dke_yp;
0653                 MMXS_t = id_MUMPSMEX_P.SOL; 
0654                 flaginv_f = flag_MUMPSMEX; 
0655             elseif dkeparam.invproc == -3,%SUPERLU 3.0 sequential MEX file (direct inversion)
0656                 [MMXS_t] = lusolve(MMXPC_P_f_t,MMXR_t);flaginv_f = 0;
0657             elseif dkeparam.invproc == -4,%Solving with PETSc
0658                 [MMXS_t,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_P_f_t,full(MMXR_t),display_mode,dkepath);
0659             end
0660         end
0661         %
0662         if flaginv_f == 2 & dkeparam.invproc > 1,
0663             error('WARNING: Matrix ill-conditionned. No inversion possible for the zero order Fokker-Planck equation. Try to reduce drop tolerance parameter for LU factorization...!');
0664         elseif flaginv_f ~= 0 & dkeparam.invproc == -1,
0665             error(['WARNING: Error ',int2str(flaginv_f),' in MUMPS matrix inversion for the zero order Fokker-Planck equation. See MUMPS documentation !']);
0666         elseif flaginv_f ~= 0 & dkeparam.invproc == -4,
0667             error(['WARNING: Error ',int2str(flaginv_f),' in PETSc matrix inversion for the zero order Fokker-Planck equation. See PETSc documentation !']);
0668         end
0669         %
0670         MMXf0_t = MMXfinit_t + MMXS_t;
0671         %
0672         for ir = 1:nr,
0673             Xf0_t = reshape(MMXf0_t(sm_f(ir)+1:sm_f(ir+1)),nmhu - ns_f(ir),npn)';
0674             %
0675             if bounce_mode == 0,
0676                 Xf0 = Xf0_t;
0677             elseif bounce_mode == 1,
0678                 Xf0 = zeros(npn,nmhu);
0679                 Xf0(:,1:nmhu/2 - ns_f(ir)) = Xf0_t(:,1:nmhu/2 - ns_f(ir));
0680                 Xf0(:,nmhu/2 + 1:nmhu) = Xf0_t(:,nmhu/2 + 1 - ns_f(ir):nmhu - ns_f(ir));
0681                 Xf0(:,nmhu/2 + 1 - ns_f(ir):nmhu/2) = fliplr(Xf0_t(:,nmhu/2 + 1 - ns_f(ir):nmhu/2));
0682             end
0683             %
0684             % compensate for avalanche source term to maintain density conservation
0685             %
0686             Xf0 = Xf0 - XXfM(:,:,ir)*dtn*xnormavalanches(ir)/xne_norm(ir);
0687             %
0688             % normf0 must be FS averaged to be compared with xne_norm
0689             %
0690             normf0(ir,it_f+1) = 2*pi*(xqtilde(ir)./xqhat(ir))*integral_dke_jd(dmhu,integral_dke_jd(dpn,XXlambda(:,:,ir).*Xf0.*Xpn2,1),1);
0691             curr0(ir,it_f+1) = 2*pi*(xq(ir)./xqbar(ir))*integral_dke_jd(dmhu,integral_dke_jd(dpn,Xf0.*Xpn2.*Xpn.*Xmhu./Xgamma,1),1);
0692             %
0693             prec_curr0(ir,it_f) = (curr0(ir,it_f+1) - curr0(ir,1))/curr0(ir,1);
0694             prec_curr(ir,it_f) = (curr0(ir,it_f+1) - curr0(ir,it_f))/curr0(ir,it_f);
0695             prec_normf0(ir,it_f) = (normf0(ir,it_f+1) - normf0(ir,1))/normf0(ir,1);
0696             prec_normf(ir,it_f) = (normf0(ir,it_f+1) - normf0(ir,it_f))/normf0(ir,it_f);
0697             %
0698             if real(dkeparam.norm_mode_f) >= 1,% use real(dkeparam.norm_mode_f) == 0 if physical losses or gains are permitted (runaway, ripple, avalanches)
0699                 %                          in that case, density conservation can be ensured in fmoments_dke_yp using boundary losses calculations
0700                 if real(dkeparam.norm_mode_f) == 1,
0701                     Xf0 = Xf0*xne_norm(ir)/normf0(ir,it_f+1);%Normalization to the maxwellian local density. Any loss or gain is assumed to be numerical.
0702                 elseif real(dkeparam.norm_mode_f) == 2,
0703                     Xf0 = Xf0*normf0(ir,1)/normf0(ir,it_f+1);%Normalization to the initial local density. Any loss or gain is assumed to be numerical.
0704                 end
0705                 %
0706                 if isfield(dkeparam,'convtest'),%Absolute residu calculation for convergence tests
0707                     residu_f(ir,it_f) = sqrt(2*pi*integral_dke_jd(dpn',pn2'.*integral_dke_jd(dmhu,(Xf0 - XXf0_ref(:,:,ir)).^2.*XXf0_ref(:,:,ir),2),1)/normf0(ir,it_f+1)/dtn/dtn);%As done by J.P. Bizarro et al., Nuclear Fusion, 37 (1997) 1509
0708                 else
0709                     residu_f(ir,it_f) = sqrt(2*pi*integral_dke_jd(dpn',pn2'.*integral_dke_jd(dmhu,(Xf0 - XXf0(:,:,ir)).^2.*XXf0(:,:,ir),2),1)/normf0(ir,it_f+1)/dtn/dtn);%As done by J.P. Bizarro et al., Nuclear Fusion, 37 (1997) 1509
0710                 end
0711                 XXf0(:,:,ir) = Xf0;
0712                 %
0713             else
0714                 residu_f(ir,it_f) = sqrt(2*pi*integral_dke_jd(dpn',pn2'.*integral_dke_jd(dmhu,(Xf0-XXf0(:,:,ir)).^2.*XXf0(:,:,ir),2),1)/normf0(ir,it_f+1)/dtn/dtn);%As done by J.P. Bizarro et al., Nuclear Fusion, 37 (1997) 1509
0715                 XXf0(:,:,ir) = Xf0;%The default value in principle. If OK, tells that the numerical conservative scheme is correct
0716             end
0717             %
0718             if it_f > 1,
0719                 prec_residu_f(ir,it_f) = abs((residu_f(ir,it_f) - residu_f(ir,it_f-1))/residu_f(ir,it_f));%For radial transport calculations
0720             end
0721             %
0722             if display_mode >= 1,
0723                 if nit_f > 1,
0724                     disp(' ');                
0725                     info_dke_yp(2,['Quasilinear iteration : ',int2str(it_rf),', Legendre iteration : ',int2str(it_f),' / ',int2str(nit_f),', Residu = ',num2str(residu_f(ir,it_f)),' / ',num2str(prec0_f)]);
0726                     info_dke_yp(2,['   Radius index number: ',int2str(ir),'; Normalized poloidal magnetic flux: ',num2str(xpsin(ir)),'; Normalized radius on LFS midplane: ',num2str(xrho(ir))]);
0727                     disp(['                - Density [it-1]: ',num2str(normf0(ir,it_f)),' - [it]: ',num2str(normf0(ir,it_f+1)),'; Variarion [t-1] : ',num2str(prec_normf0(ir,it_f)),'; Variarion [it-1] : ',num2str(prec_normf(ir,it_f))]);               
0728                     disp(['                - Current [it-1]: ',num2str(curr0(ir,it_f)),' - [it]: ',num2str(curr0(ir,it_f+1)),'; Variarion [t-1] : ',num2str(prec_curr0(ir,it_f)),' Variarion [it-1] : ',num2str(prec_curr(ir,it_f))]);
0729                 else
0730                     disp(' ');                
0731                     info_dke_yp(2,['Quasilinear iteration : ',int2str(it_rf)]);
0732                     info_dke_yp(2,['   Radius index number: ',int2str(ir),'; Normalized poloidal magnetic flux: ',num2str(xpsin(ir)),'; Normalized radius on LFS midplane: ',num2str(xrho(ir))]);
0733                     disp(['                - Density [it-1]: ',num2str(normf0(ir,it_f)),' - [it]: ',num2str(normf0(ir,it_f+1)),'; Variarion [t-1] : ',num2str(prec_normf0(ir,it_f))]);               
0734                     disp(['                - Current [it-1]: ',num2str(curr0(ir,it_f)),' - [it]: ',num2str(curr0(ir,it_f+1)),'; Variarion [t-1] : ',num2str(prec_curr0(ir,it_f))]);
0735                 end
0736             end
0737             %            %
0738             if display_mode == 2,
0739                 if dke_mode == 1,
0740                     if ir == ir_display-1,
0741                         figure(Fig_FPEconvstatus),subplot(321),semilogy(-pn,XXfM(:,2,ir),'g--',pn,XXfM(:,nmhu-1,ir),'g--',-pn,Xf1(:,2),'r-',pn,Xf1(:,nmhu-1),'r-');%f(perp = 0)
0742                         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
0743                         title(['Iteration: ',int2str(it_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
0744                         hold on;zoom on
0745                         figure(Fig_FPEconvstatus),subplot(322),semilogy([1:it_f],residu_f(ir,1:it_f),'r+',[1:it_f],max(residu_f(:,1:it_f)),'gx',[1:it_f],prec0_f*ones(size([1:it_f])),'b--');
0746                         axis('square');axis([0,it_f,prec0_f/10,1e-4]);
0747                         title(['Convergence - f (rho=',num2str(xrho(ir)),')']);xlabel('Iteration');ylabel('Residu')
0748                         hold on;zoom on
0749                         drawnow;
0750                     elseif ir == ir_display,
0751                         figure(Fig_FPEconvstatus),subplot(323),semilogy(-pn,XXfM(:,2,ir),'g--',pn,XXfM(:,nmhu-1,ir),'g--',-pn,Xf1(:,2),'r-',pn,Xf1(:,nmhu-1),'r-');%f(perp = 0)
0752                         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
0753                         title(['Iteration: ',int2str(it_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
0754                         hold on;zoom on
0755                         figure(Fig_FPEconvstatus),subplot(324),semilogy([1:it_f],residu_f(ir,1:it_f),'r+',[1:it_f],max(residu_f(:,1:it_f)),'gx',[1:it_f],prec0_f*ones(size([1:it_f])),'b--');
0756                         axis('square');axis([0,it_f,prec0_f/10,1e-4]);
0757                         title(['Convergence - f (rho=',num2str(xrho(ir)),')']);xlabel('Iteration');ylabel('Residu')
0758                         hold on;zoom on
0759                         drawnow;
0760                     elseif ir == ir_display+1,
0761                         figure(Fig_FPEconvstatus),subplot(325),semilogy(-pn,XXfM(:,2,ir),'g--',pn,XXfM(:,nmhu-1,ir),'g--',-pn,Xf1(:,2),'r-',pn,Xf1(:,nmhu-1),'r-');%f(perp = 0)
0762                         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
0763                         title(['Iteration: ',int2str(it_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
0764                         hold on;zoom on
0765                         figure(Fig_FPEconvstatus),subplot(326),semilogy([1:it_f],residu_f(ir,1:it_f),'r+',[1:it_f],max(residu_f(:,1:it_f)),'gx',[1:it_f],prec0_f*ones(size([1:it_f])),'b--');
0766                         axis('square');axis([0,it_f,prec0_f/10,1e-4]);
0767                         title(['Convergence - f (rho=',num2str(xrho(ir)),')']);xlabel('Iteration');ylabel('Residu')
0768                         hold on;zoom on
0769                         drawnow;
0770                     end
0771                 else
0772                     if ir == ir_display,
0773                         figure(Fig_FPEconvstatus),subplot(121),semilogy(-pn,XXfM(:,2,ir),'g--',pn,XXfM(:,nmhu-1,ir),'g--',-pn,Xf0(:,2),'r-',pn,Xf0(:,nmhu-1),'r-');%f(perp = 0)
0774                         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
0775                         title(['Iteration: ',int2str(it_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
0776                         hold on;
0777                         figure(Fig_FPEconvstatus),subplot(122),semilogy([1:it_f],residu_f(ir,1:it_f),'r+',[1:it_f],max(residu_f(:,1:it_f)),'gx',[1:it_f],prec0_f*ones(size([1:it_f])),'b--');
0778                         axis('square');axis([0,it_f,prec0_f/10,1e-4]);
0779                         title(['Convergence - f (rho=',num2str(xrho(ir)),')']);xlabel('Iteration');ylabel('Residu')
0780                         hold on;
0781                         drawnow;
0782                     end
0783                 end
0784             end
0785         end
0786         %
0787         %
0788         if display_mode == -1,    
0789             %
0790             figure(1),plot(xrho,curr0(:,it_f+1));hold on
0791             figure(2),semilogy(it_f+1,max(residu_f(:,it_f)),'r+');hold on
0792             %
0793             [ir_max] =  find(residu_f(:,it_f) == max(residu_f(:,it_f)));
0794             disp(' ');                
0795             info_dke_yp(2,['Quasilinear iteration : ',int2str(it_rf),', Fokker-Planck iteration : ',int2str(it_f),' / ',int2str(nit_f),', Residu max.= ',num2str(residu_f(ir_max,it_f)),' / ',num2str(prec0_f)]);
0796             info_dke_yp(2,['   Normalized poloidal magnetic flux: ',num2str(xpsin(ir_max))]);
0797             info_dke_yp(2,['   Normalized radius on LFS midplane: ',num2str(xrho(ir_max)),' / Radius index number: ',int2str(ir_max)]);
0798             disp(['                - Density [t-1]: ',num2str(normf0(ir_max,it_f)),' - [t]: ',num2str(normf0(ir_max,it_f+1)),' / Absolute variation [t-t0]: ',num2str(prec_normf0(ir_max,it_f)),' / Relative variation [dt]: ',num2str(prec_normf(ir_max,it_f))]);               
0799             disp(['                - Current [t-1]: ',num2str(curr0(ir_max,it_f)),' - [t]: ',num2str(curr0(ir_max,it_f+1)),' / Current relative accuracy [dt]: ',num2str(prec_curr(ir_max,it_f))]);
0800         end
0801         %
0802      %   if abs(max(residu_f(:,it_f))) > 1e-6,
0803      %       info_dke_yp(2,['WARNING: Divergence in the calculation of f0 after ',int2str(it_f),' iterations. Time step likely too large !']);
0804      %       break;
0805      %   end
0806         %
0807         if all(residu_f(:,it_f) <= prec0_f),
0808             if sum(prec_residu_f(:,it_f)) < prec0_f,info_dke_yp(2,['WARNING: stationnary residu but convergence not achieved for f0 after ',int2str(it_f),' iterations. Time step likely too large !']);end
0809             if abs(display_mode) > 0,
0810                 info_dke_yp(2,['Successfull convergence achieved for f0 after ',int2str(it_f),' iterations !']);
0811             end
0812             break;
0813         end
0814     end
0815     %
0816     if it_f == nit_f && prec0_f >= 0,info_dke_yp(2,['WARNING: convergence not achieved for f0 after ',int2str(nit_f),' iterations. Time step too small !']);end
0817     %
0818 end
0819 %
0820 if dkeparam.invproc == -2,%MUMPS FORTAN 95 parallel processing package (inversion in MEX file, with LU matrix preconditionning)
0821    job_MUMPSMEX = -2;MUMPSMEX_dke_yp;%Clear MUMPSMEX memory
0822 end
0823 if abs(display_mode) > 0,
0824     info_dke_yp(2,['Convergence achieved for zero order Fokker-Planck calculations in ',num2str(etime(clock,time0)),' (s)']);
0825 end
0826 if dkeparam.invproc == -1 && abs(display_mode) > 0,info_dke_yp(2,['Number of processors used for MUMPS matrix inversion : ',int2str(dkeparam.MUMPSparam.nproc)]);end%MUMPS FORTAN 95 parallel processing package
0827 %
0828 clear ir fM XfM MMXPC_f_t MMXPDI_f_t Xmask_f0 MMXL_f_t MMXU_f_t
0829 clear Xf0 Xf0_t Xq Xq_t Xf1 Xf1_t Xf1_tp Xf1_tm Xf1_ts XR_t XS_t XY_t
0830 clear f1 xf1 I1 I2 XI 
0831 %
0832 if display_mode == 2,
0833     figure(Fig_FPEconvstatus),subplot(121),hold off;
0834     figure(Fig_FPEconvstatus),subplot(122),hold off;
0835 end
0836 %
0837 time0 = clock;
0838 %
0839 if dke_mode == 1,
0840     if display_mode == 2,
0841         Fig_DKEconvstatus = figure('Name','Drift kinetic equation convergence status');
0842     end    
0843 %
0844     MMXf0_tp = [];
0845     MMXf0 = [];
0846     MMXI_tp = [];
0847     MMXq_tp_t = []; 
0848 %
0849 %Calculation of the right handside term of the first order Fokker-Planck equation
0850 %
0851     for ir_dke = 1:nr_dke,
0852         jr1 = rdke(ir_dke)-1;
0853         jr2 = rdke(ir_dke);
0854         jr3 = rdke(ir_dke)+1;
0855 %
0856         Xf0m = XXf0(:,:,jr1);%Lower radial value for gradient calculation
0857         Xf0  = XXf0(:,:,jr2);%Central radial value for gradient calculation (where ftp and g are calculated)
0858         Xf0p = XXf0(:,:,jr3);%Upper radial value for gradient calculation
0859 %
0860         XfMm = XXfM(:,:,jr1);%Lower radial value for gradient calculation
0861         XfM = XXfM(:,:,jr2);%Central radial value for gradient calculation (where ftp and g are calculated)
0862         XfMp = XXfM(:,:,jr3);%Upper radial value for gradient calculation
0863 %
0864         dpsim = (xpsin(jr1) - xpsin(jr2));%Radial step for parabolic gradient calculations (dke)
0865         dpsi = (xpsin(jr3) - xpsin(jr1));%Radial step for parabolic gradient calculations (dke)
0866         dpsip = (xpsin(jr3) - xpsin(jr2));%Radial step for parabolic gradient calculations (dke)
0867 %
0868         %if dkeparam.coll_mode == 3,%Normalized gradients of ne and Te for the Lorentz model (test mode only, see documentation)
0869             xdlnnedpsi(jr2) = (-dpsim*xne_norm(jr3)/dpsi/dpsip - (dpsip + dpsim)*xne_norm(jr2)/dpsip/dpsim + dpsip*xne_norm(jr1)/dpsi/dpsim)/xne_norm(jr2);
0870             xdlnTedpsi(jr2) = (-dpsim*xTe_norm(jr3)/dpsi/dpsip - (dpsip + dpsim)*xTe_norm(jr2)/dpsip/dpsim + dpsip*xTe_norm(jr1)/dpsi/dpsim)/xTe_norm(jr2);
0871         %end
0872 %
0873         Xdf0dpsi = -dpsim*Xf0p/dpsi/dpsip - (dpsip + dpsim)*Xf0/dpsip/dpsim + dpsip*Xf0m/dpsi/dpsim;%Calculation of the radial gradient of f0 (parabolic interpolation)
0874         XXf0_tp(:,:,jr2) = -xftp_norm(jr2)*(Xpn.*Xmhu.*Xdf0dpsi);
0875 %
0876 % Initial guess for the iterative loop (Lorentz model)
0877 %
0878         XdfMdpsi = -dpsim*XfMp/dpsi/dpsip - (dpsip + dpsim)*XfM/dpsip/dpsim + dpsip*XfMm/dpsi/dpsim;%Calculation of the radial gradient of fM (parabolic interpolation)
0879         XXf0_g_th(:,:,jr2) = +xftp_norm(jr2)*(xdlnnedpsi(jr2) + xdlnTedpsi(jr2)*(Xpn2/xTe_norm(jr2)/2 - 1.5)).*Xpn.*(Xmhu./abs(Xmhu)).*XXfM(:,:,jr2).*abs(XXILor(:,:,jr2));%
0880 %
0881 %First order Legendre correction for ftp
0882 %
0883         XI_tp = zeros(npn,nmhu);
0884         if dkeparam.coll_mode == 2,%Belaiev-Budker collision operator
0885             f1_tp = 1.5*integral_dke_jd(dmhu,XXf0_tp(:,:,jr2).*XXlambda_b_p2(:,:,jr2).*(ones(npn,1)*mhu),2)';
0886             [I1_tp,I2_tp] = firstordercollop_dke_yp(f1_tp,betath_ref,...
0887                      masku,maskl,...
0888                      dpn,xpn,xpn2,xpn3,xsigma,xgamma,xz,xz2,pn,pn2,pn3,sigma,gamma,z,z2,...
0889                      xJ1,xJ2,xJ3,xJ4,J1,J2,J3,J4,...
0890                      xTe_norm(jr2));
0891             XI_tp = 4*pi*((XXfM(:,1,jr2)'.*(f1_tp + I1_tp./pn2 + I2_tp.*pn)./gamma)'*mhu).*XXheaviside(:,:,jr2)./XXlambda(:,:,jr2);
0892         else
0893             XI_tp = zeros(npn,nmhu);
0894         end
0895 %
0896         MMXf0_tp = [MMXf0_tp;reshape(XXf0_tp(:,:,jr2)',nmhu*npn,1)];        
0897         MMXI_tp = [MMXI_tp;reshape(XI_tp'.*Xpn2',nmhu*npn,1)];%Multiplication by the Jacobian of momentum coordinates
0898     end
0899 %
0900     MMXq_tp = -MMXP_tp*MMXf0_tp + MMXI_tp;%****** WARNING : en mode Lorentz et en multipliant MMXP_tp_t*MMXf0_tp par 0.5 on retrouve le bon niveau de courant ****
0901 %    MMXq_tp_c = -MMH_tp*reshape(permute(XXf0,[2,1,3]),npn*nmhu*nr,1,1);
0902 %
0903     for ir_dke = 1:nr_dke,%Remove trapped region at jr2
0904         Xq_tp = reshape(MMXq_tp(sm_tp(ir_dke)+1:sm_tp(ir_dke+1)),nmhu,npn)';
0905 %        Xq_tp_c = reshape(MMXq_tp_c(npn*nmhu*rdke(ir_dke)+1:npn*nmhu*(rdke(ir_dke)+1)),nmhu,npn)';
0906         if boundary_mode_g ~= 0,
0907             Xq_tp(1:boundary_mode_g,:) = zeros(1:boundary_mode_g,nmhu);%Boundary conditions at p = 0 enforced for g
0908 %            Xq_tp_c(1:boundary_mode_g,:) = zeros(1:boundary_mode_g,nmhu);%Boundary conditions at p = 0 enforced for g
0909         end
0910  %      MMXq_tp_t = [MMXq_tp_t;reshape(Xq_tp(:,Zmask_g{ir_dke})' + Xq_tp_c(:,Zmask_g{ir_dke})',(nmhu - ns_g(ir_dke))*npn,1)];%Matrix reduction since the trapped region is removed
0911         MMXq_tp_t = [MMXq_tp_t;reshape(Xq_tp(:,Zmask_g{ir_dke})',(nmhu - ns_g(ir_dke))*npn,1)];%Matrix reduction since the trapped region is removed
0912     end   
0913 %
0914 %Resolution of the first order Fokker-Planck equation
0915 %
0916     XXf0_g = XXf0_g_th;%Initial guess (thermal limit, large aspect ratio, Lorentz model)
0917     for it_g = 2:nit_g, 
0918         if sum(residu_g(:,it_g-1) >= prec0_g) >= 1 & sum(prec_residu_g(:,it_g-1)) >= prec0_g | it_g < 5,
0919 %
0920             MMXf0_g_t = [];
0921             MMXq_g_t = [];
0922             for ir_dke = 1:nr_dke,
0923                 jr2 = rdke(ir_dke);
0924 %
0925                 Xq_g_t = XXI_g(:,Zmask_g{ir_dke},jr2).*Xpn2(:,Zmask_g{ir_dke});%First order expansion of the Beliaev-Budker collision operator (current drive problem only)
0926                 if boundary_mode_g ~= 0,
0927                         Xq_g_t(1:boundary_mode_g,:) = XXf0_g_th(1:boundary_mode_g,Zmask_g{ir_dke},jr2);%Boundary conditions at p = 0 (large aspect ratio limit, Lorentz model)
0928                 end
0929 %
0930                 MMXf0_g_t = [MMXf0_g_t;reshape(XXf0_g(:,Zmask_g{ir_dke},jr2)',(nmhu - ns_g(ir_dke))*npn,1)];
0931                 MMXq_g_t = [MMXq_g_t;reshape(Xq_g_t',(nmhu - ns_g(ir_dke))*npn,1)];
0932             end
0933 %
0934             MMXR_t = sparse(MMXPDI_g_t*(MMXq_g_t + MMXq_tp_t) - MMXPC_g_t*MMXf0_g_t);%Force sparse tyep for the RHS, optmized for Matlab
0935 %
0936             if dkeparam.invproc == 0,%Matlab direct inversion
0937                MMXS_t = MMXR_t\MMXPC_g_t;
0938             elseif dkeparam.invproc == 1,%cgs iterative method
0939                [MMXS_t,flaginv_g] = cgs(MMXPC_g_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_g_t,MMXU_g_t);
0940             elseif dkeparam.invproc == 2,%bicg iterative method
0941                [MMXS_t,flaginv_g] = bicg(MMXPC_g_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_g_t,MMXU_g_t); 
0942             elseif dkeparam.invproc == 3,%qmr iterative method
0943                [MMXS_t,flaginv_g] = qmr(MMXPC_g_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_g_t,MMXU_g_t);
0944             elseif dkeparam.invproc == -1,%MUMPS FORTAN 95 parallel processing package (direct inversion)
0945                [MMXS_t,flaginv_g] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_g_t,full(MMXR_t),display_mode,dkepath);
0946             elseif dkeparam.invproc == -4,%Solving with PETSc
0947                [MMXS_t,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),display_mode,dkepath);
0948             end
0949 %
0950             if flaginv_g == 2 & dkeparam.invproc > 1,
0951                error('WARNING: Matrix ill-conditionned. No inversion possible for the first order Fokker-Planck equation. Try to reduce drop tolerance parameter for LU factorization...!');
0952             elseif flaginv_g ~= 0 & dkeparam.invproc == -1,
0953                error(['WARNING: Error ',int2str(flaginv_g),' in MUMPS matrix inversion for the first order Fokker-Planck equation. See MUMPS documentation !']);
0954             elseif flaginv_g ~= 0 & dkeparam.invproc == -4,
0955                error(['WARNING: Error ',int2str(flaginv_g),' in PETSc matrix inversion for the first order Fokker-Planck equation. See PETSc documentation !']);
0956             end
0957 %
0958             MMXf1_g_t = MMXf0_g_t + MMXS_t;
0959 %
0960             for ir_dke = 1:nr_dke,
0961                 jr2 = rdke(ir_dke);
0962                 Xf1_g_t = reshape(MMXf1_g_t(sm_g(ir_dke)+1:sm_g(ir_dke+1)),nmhu - ns_g(ir_dke),npn)';%Split vector into sub-matrices
0963 %
0964                 Xf1_g = zeros(npn,nmhu);
0965                 Xf1_g(1:npn,1:(nmhu - ns_g(ir_dke))/2) = Xf1_g_t(1:npn,1:(nmhu - ns_g(ir_dke))/2);%Expansion of the function g on the pitch-angle grid
0966                 Xf1_g(1:npn,(nmhu + ns_g(ir_dke))/2 + 1:nmhu) = Xf1_g_t(1:npn,(nmhu - ns_g(ir_dke))/2 + 1:nmhu - ns_g(ir_dke));%Expansion of the function g on the pitch-angle grid
0967 %
0968                 norm1_g(ir_dke,it_g) = 2*pi*integral_dke_jd(dpn',pn2'.*integral_dke_jd(dmhu,Xf1_g,2),1);
0969                 residu_g(ir_dke,it_g) = sqrt(abs(2*pi*integral_dke_jd(dpn',pn2'.*integral_dke_jd(dmhu,(Xf1_g - XXf0_g(:,:,jr2)).^2.*XXf0_g(:,:,jr2),2),1)/norm1_g(ir_dke,it_g)));
0970                 prec_residu_g(ir_dke,it_g) = abs((residu_g(ir_dke,it_g) - residu_g(ir_dke,it_g-1))/residu_g(ir_dke,it_g));              
0971 %
0972                 if display_mode >= 1,   
0973                     info_dke_yp(2,['Drift kinetic iteration : ',int2str(it_g),' / ',int2str(nit_g),', Residu = ',num2str(residu_g(ir_dke,it_g))]);
0974                     info_dke_yp(2,['   Normalized poloidal magnetic flux: ',num2str(xpsin(jr2))]);
0975                     info_dke_yp(2,['      Normalized radius on LFS midplane: ',num2str(xrho(jr2)),' / Radius index number: ',int2str(jr2)]);
0976                 end
0977 %
0978                 if jr2 == ir_display & display_mode == 2,
0979                     figure(Fig_DKEconvstatus),subplot(121),plot(-pn,Xf1_g(:,1),'r-',pn,Xf1_g(:,nmhu),'r-');%g(perp = 0)
0980                     axis('square');set(gca,'xlim',[-10,10]);
0981                     title(['Iteration: ',int2str(it_g)]);xlabel('pn-ref');ylabel('g(pperp=0)')
0982                     hold on;zoom on
0983                     figure(Fig_DKEconvstatus),subplot(122),semilogy([1:it_g],residu_g(ir_dke,1:it_g),'r+',[1:it_g],max(residu_g(:,1:it_g)),'gx',[1:it_g],prec0_g*ones(size([1:it_g])),'b--');
0984                     axis('square');axis([0,it_g,prec0_g/10,1e-4]);
0985                     title(['Convergence - g (rho=',num2str(xrho(jr2)),')']);xlabel('Iteration');ylabel('Residu')
0986                     hold on;zoom on
0987                     drawnow;
0988                 end
0989 %
0990 %First order Legendre correction for g
0991 %
0992                 if dkeparam.coll_mode == 2,%Belaiev-Budler collision operator
0993                     f1_g = 1.5*integral_dke_jd(dmhu,Xf1_g.*(ones(npn,1)*mhu),2)';
0994                     [I1_g,I2_g] = firstordercollop_dke_yp(f1_g,betath_ref,...
0995                              masku,maskl,...
0996                              dpn,xpn,xpn2,xpn3,xsigma,xgamma,xz,xz2,pn,pn2,pn3,sigma,gamma,z,z2,...
0997                              xJ1,xJ2,xJ3,xJ4,J1,J2,J3,J4,...
0998                              xTe_norm(jr2));
0999                          XI_g = 4*pi*((XXfM(:,1,jr2)'.*(f1_g + I1_g./pn2 + I2_g.*pn)./gamma)'*mhu).*XXRlambda_b_p1p1(:,:,jr2);   
1000                          XXI_g(:,:,jr2) = XI_g;
1001                 else
1002                     XXI_g(:,:,jr2) = zeros(npn,nmhu);
1003                 end
1004 %
1005                 XXf0_g(:,:,jr2) = Xf1_g;
1006             end
1007         else
1008             if sum(prec_residu_g(:,it_g-1)) < prec0_g,info_dke_yp(2,['WARNING: stationnary residu but convergence not achieved for g after ',int2str(it_g),' iterations. Time step likely too large !']);end
1009             if sum(residu_g(:,it_g-1) >= prec0_g) == 0 & display_mode > 0,
1010                 info_dke_yp(2,['Successfull convergence achieved for g after ',int2str(it_g),' iterations !']);
1011             end
1012             break;  
1013         end
1014         if it_g == nit_g,info_dke_yp(2,['WARNING: convergence not achieved for g after ',int2str(nit_g),' iterations.']);end
1015     end
1016 %
1017     if display_mode >= 1,info_dke_yp(2,['Convergence achieved for first order Fokker-Planck calculations in ',num2str(etime(clock,time0)),' (s)']);end
1018     if dkeparam.invproc == -1 & display_mode >= 1,info_dke_yp(2,['Number of processors used for MUMPS matrix inversion : ',int2str(dkeparam.MUMPSparam.nproc)]);end%MUMPS FORTAN 95 parallel processing package
1019 else
1020     XXf0_tp = zeros(size(XXf0));
1021     XXf0_g = zeros(size(XXf0));
1022 end
1023 %
1024 if dkeparam.coll_mode == 3 & dke_mode == 1,%Lorentz model (test mode only, see documentation)
1025     for ir_dke = 1:nr_dke,
1026         jr2 = rdke(ir_dke);
1027         XXf0_L_tp(:,:,jr2) = -xftp_norm(jr2)*(xdlnnedpsi(jr2) + xdlnTedpsi(jr2)*(Xpn2/xTe_norm(jr2)/2 - 1.5)).*Xpn.*Xmhu.*XXfM(:,:,jr2);
1028         XXf0_L_g(:,:,jr2) =  xftp_norm(jr2)*(xdlnnedpsi(jr2) + xdlnTedpsi(jr2)*(Xpn2/xTe_norm(jr2)/2 - 1.5)).*Xpn.*(Xmhu./abs(Xmhu)).*XXfM(:,:,jr2).*abs(XXILor(:,:,jr2));%Signe ********
1029     end 
1030 end 
1031 %
1032 if display_mode == 2,
1033     if dke_mode == 1,
1034 %
1035         XfM = XXfM(:,:,ir_display-1);
1036         Xf0 = XXf0(:,:,ir_display-1);
1037         Fig_Finalconvstatus = figure('Name','Final convergence status');subplot(321),semilogy(-pn,XfM(:,2),'g--',pn,XfM(:,nmhu-1),'g--',-pn,Xf0(:,2),'r-',pn,Xf0(:,nmhu-1),'r-');%f(perp = 0)
1038         axis('square');axis([-max(pn),max(pn),1e-20,10]);zoom on
1039         title(['Iteration: ',int2str(nit_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
1040         hold on;zoom on
1041         figure(Fig_Finalconvstatus),subplot(322),plot([1:nit_f],curr1(ir_display-1,:),'r-');
1042         axis('square');if max(max(curr1(ir_display-1,:))) ~= 0,axis([0,nit_f,0,1.1*max(max(curr1(ir_display-1,:)))]);end
1043         title('Convergence');xlabel('Iteration');ylabel('Current')
1044         hold on;zoom on
1045         drawnow;
1046 %
1047         XfM = XXfM(:,:,ir_display);
1048         Xf0 = XXf0(:,:,ir_display);
1049         figure(Fig_Finalconvstatus),subplot(323),semilogy(-pn,XfM(:,2),'g--',pn,XfM(:,nmhu-1),'g--',-pn,Xf0(:,2),'r-',pn,Xf0(:,nmhu-1),'r-');%f(perp = 0)
1050         axis('square');axis([-max(pn),max(pn),1e-20,10]);zoom on
1051         title(['Iteration: ',int2str(nit_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
1052         hold on;zoom on
1053         figure(Fig_Finalconvstatus),subplot(324),plot([1:nit_f],curr1(ir_display,:),'r-');
1054         axis('square');if max(max(curr1(ir_display,:))) ~= 0,axis([0,nit_f,0,1.1*max(max(curr1(ir_display,:)))]);end
1055         title('Convergence');xlabel('Iteration');ylabel('Current')
1056         hold on;zoom on
1057         drawnow;
1058 %
1059         XfM = XXfM(:,:,ir_display+1);
1060         Xf0 = XXf0(:,:,ir_display+1);
1061         figure(Fig_Finalconvstatus),subplot(325),semilogy(-pn,XfM(:,2),'g--',pn,XfM(:,nmhu-1),'g--',-pn,Xf0(:,2),'r-',pn,Xf0(:,nmhu-1),'r-');%f(perp = 0)
1062         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
1063         title(['Iteration: ',int2str(nit_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
1064         hold on;zoom on
1065         figure(Fig_Finalconvstatus),subplot(326),plot([1:nit_f],curr1(ir_display+1,:),'r-');
1066         axis('square');if max(max(curr1(ir_display+1,:))) ~= 0,axis([0,nit_f,0,1.1*max(max(curr1(ir_display+1,:)))]);end
1067         title('Convergence');xlabel('Iteration');ylabel('Current')
1068         hold on;zoom on
1069         drawnow;
1070     else
1071         XfM = XXfM(:,:,ir_display);
1072         Xf0 = XXf0(:,:,ir_display);
1073 %
1074         Fig_Finalconvstatus = figure('Name','Final convergence status');
1075         figure(Fig_Finalconvstatus),subplot(121),semilogy(-pn,XfM(:,2),'g--',pn,XfM(:,nmhu-1),'g--',-pn,Xf0(:,2),'r-',pn,Xf0(:,nmhu-1),'r-');%f(perp = 0)
1076         axis('square');axis([-max(pn),max(pn),1e-20,10]),zoom on
1077         title(['Iteration: ',int2str(nit_f)]);xlabel('pn-ref');ylabel('f(pperp=0)')
1078         hold on;
1079         figure(Fig_Finalconvstatus),subplot(122),plot([1:nit_f],curr1(ir_display,:),'r-');
1080         axis('square');if max(max(curr1(ir_display,:))) ~= 0,axis([0,nit_f,0,1.1*max(max(curr1(ir_display,:)))]);end
1081         title('Current convergence');xlabel('Iteration');ylabel('Current')
1082         hold on;
1083         drawnow;
1084     end
1085 end
1086 %
1087 time_DKE = etime(clock,time00);
1088 %
1089 if display_mode >= 1,info_dke_yp(2,['Full DKE equation solved in ',num2str(time_DKE),' (s)']);end