rundke

Script for running the DKE solver (can be modified by the user for specific simulations)
by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)

0001 %Script for running the DKE solver (can be modified by the user for specific simulations)
0002 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)
0003 %
0004 clear all
0005 clear mex
0006 clear functions
0007 close all
0008 warning off
0009 global nfig
0010 %
0011 p_opt = 2;
0012 %
0013 permission = test_permissions_yp;
0014 %
0015 if ~permission 
0016     disp('Please move the script to a local folder where you have write permission before to run it')
0017     return;
0018 end
0019 %
0020 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0021 %
0022 id_simul = 'LH_karney_num';%Simulation ID
0023 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0024 %
0025 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0026 rho_S = [0.5];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0027 %
0028 id_path = '';%For all paths used by DKE solver
0029 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0030 %
0031 id_equil = 'TScyl';%For plasma equilibrium
0032 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0033 %
0034 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0035 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0036 %
0037 id_display = 'NO_DISPLAY';%For output code display
0038 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0039 %
0040 id_ohm = '';%For Ohmic electric contribution
0041 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0042 %
0043 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0044 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0045 %
0046 id_transpfaste = '';%For fast electron radial transport
0047 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0048 %
0049 id_ripple = '';%For fast electron magnetic ripple losses
0050 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0051 %
0052 %************************************************************************************************************************************
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %
0056 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple);
0057 %
0058 %************************************************************************************************************************************
0059 %
0060 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0061 %Option parameter for cross-comparison between old LH code:
0062 %    - (1): 1/vpar dependence
0063 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0064 wavestruct.opt_lh = 2; % [1,1]
0065 %
0066 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0067 % or (Npar_lh,dNpar_lh) for Gaussian shape
0068 %
0069 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0070 %
0071 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0072 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0073 %
0074 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0075 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0076 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0077 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0078 %
0079 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0080 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0081 %
0082 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0083 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0084 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0085 %
0086 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0087 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0088 %
0089 %************************************************************************************************************************************
0090 %
0091 dkeparam.boundary_mode_f = 0;%Number of points where the Maxwellian distribution is enforced from p = 0 (p=0, free conservative mode but param_inv(1) must be less than 1e-4, otherwise 1e-3 is OK most of the time. Sensitive to the number of points in p)
0092 %
0093 dkeparam.psin_S = psin_S;
0094 dkeparam.rho_S = rho_S;
0095 %
0096 % with norm_mode_f = 0, uncomment n-1 plots
0097 %
0098 dkeparam.norm_mode_f = 1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0099 %
0100 % with coll_mode = 2, uncomment nit plots
0101 %
0102 dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0103 %
0104 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0105 %
0106 dkeparam.mlmaxitinv = 20;%Maximum number of iterations for iterative built-in MatLab matrix inversion method (between 1 and 20 usualy, if empty, mlmaxitinv = min(matrixdim,20). See matlab documentation.
0107 dkeparam.mlprecinv = 1e-6;%Accuracy for the MatLab build-in matrix inversion. (default : 1e-6)
0108 %
0109 % Direct MATLAB inversion
0110 %
0111 dkeparam.invproc = 0;
0112 %
0113 [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0114 %
0115 memory_ML = dke_out.memoryLU_f_out;
0116 time_ML = dke_out.time_DKE;
0117 %
0118 nit_ML = length(dke_out.residu_f{end});
0119 %
0120 norm_ML = Znorm.x_0;
0121 curr_ML = Zcurr.x_0;
0122 Prf_ML = ZP0.x_rf_fsav;
0123 Pcoll_ML = ZP0.x_c_fsav;
0124 %
0125 disp(['Direct MatLab inversion done'])
0126 %
0127 % Exact LU factorization
0128 %
0129 dkeparam.invproc = 1;%MATLAB-CGS
0130 dkeparam.ludroptol = 0;%Full LU factorization for preconditionning
0131 %
0132 [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0133 %
0134 memory_LU = dke_out.memoryLU_f_out;
0135 time_LU = dke_out.time_DKE;
0136 %
0137 nit_LU = length(dke_out.residu_f{end});
0138 %
0139 norm_LU = Znorm.x_0;
0140 curr_LU = Zcurr.x_0;
0141 Prf_LU = ZP0.x_rf_fsav;
0142 Pcoll_LU = ZP0.x_c_fsav;
0143 %
0144 disp('Full LU calculation done')
0145 %
0146 % Incomplete LU factorization
0147 %
0148 ludroptol_list = logspace(-8,-2,7);
0149 mlprecinv_list = logspace(-8,-2,7);
0150 iludrop = 1;
0151 imlprecinv = 1;
0152 %
0153 nludroptol = length(ludroptol_list);
0154 nmlprecinv = length(mlprecinv_list);
0155 %
0156 memory_LUi_1 = NaN(nludroptol,3);
0157 time_LUi_1 = NaN(nludroptol,3);
0158 nit_LUi_1 = NaN(nludroptol,3);
0159 norm_LUi_1 = NaN(nludroptol,3);
0160 curr_LUi_1 = NaN(nludroptol,3);
0161 Prf_LUi_1 = NaN(nludroptol,3);
0162 Pcoll_LUi_1 = NaN(nludroptol,3);
0163 %
0164 memory_LUi_2 = NaN(nmlprecinv,3);
0165 time_LUi_2 = NaN(nmlprecinv,3);
0166 nit_LUi_2 = NaN(nmlprecinv,3);
0167 norm_LUi_2 = NaN(nmlprecinv,3);
0168 curr_LUi_2 = NaN(nmlprecinv,3);
0169 Prf_LUi_2 = NaN(nmlprecinv,3);
0170 Pcoll_LUi_2 = NaN(nmlprecinv,3);
0171 %
0172 for imeth = 1:3,%Numerical method (CGS, BICG, QMR)
0173     %
0174     dkeparam.invproc  = imeth;
0175     %
0176     dkeparam.mlprecinv = 1e-6;%Accuracy for the MatLab build-in matrix inversion. (default : 1e-6)
0177     %
0178     for iludroptol = 1:nludroptol,%droptolerance scan
0179         %
0180         dkeparam.ludroptol = ludroptol_list(iludroptol);
0181         %
0182         try
0183             %
0184             [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0185             %
0186             memory_LUi_1(iludroptol,imeth) = dke_out.memoryLU_f_out;
0187             time_LUi_1(iludroptol,imeth) = dke_out.time_DKE;
0188             %
0189             nit_LUi_1(iludroptol,imeth) = length(dke_out.residu_f{end});
0190             %
0191             norm_LUi_1(iludroptol,imeth) = Znorm.x_0;
0192             curr_LUi_1(iludroptol,imeth) = Zcurr.x_0;
0193             Prf_LUi_1(iludroptol,imeth) = ZP0.x_rf_fsav;
0194             Pcoll_LUi_1(iludroptol,imeth) = ZP0.x_c_fsav;
0195             %
0196             disp(['Calculation with MatLab incomplete LU calculation #',int2str(iludroptol),' with matrix inversion method #',int2str(imeth),'  done'])
0197         catch
0198             disp(['Calculation with MatLab incomplete LU calculation #',int2str(iludroptol),' with matrix inversion method #',int2str(imeth),'  failed'])
0199         end
0200         %
0201     end
0202     %
0203     dkeparam.ludroptol = 1e-6;%Accuracy for the MatLab build-in matrix inversion. (default : 1e-6)
0204     %
0205     for imlprecinv = 1:nmlprecinv,%droptolerance scan
0206         %
0207         dkeparam.mlprecinv = mlprecinv_list(imlprecinv);
0208         %
0209         try
0210             %
0211             [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0212             %
0213             memory_LUi_2(imlprecinv,imeth) = dke_out.memoryLU_f_out;
0214             time_LUi_2(imlprecinv,imeth) = dke_out.time_DKE;
0215             %
0216             nit_LUi_2(imlprecinv,imeth) = length(dke_out.residu_f{end});
0217             %
0218             norm_LUi_2(imlprecinv,imeth) = Znorm.x_0;
0219             curr_LUi_2(imlprecinv,imeth) = Zcurr.x_0;
0220             Prf_LUi_2(imlprecinv,imeth) = ZP0.x_rf_fsav;
0221             Pcoll_LUi_2(imlprecinv,imeth) = ZP0.x_c_fsav;
0222             %
0223             disp(['Calculation with MatLab incomplete LU calculation #',int2str(imlprecinv),' with matrix inversion method #',int2str(imeth),'  done'])
0224         catch
0225             disp(['Calculation with MatLab incomplete LU calculation #',int2str(imlprecinv),' with matrix inversion method #',int2str(imeth),'  failed'])
0226         end
0227         %
0228     end
0229     %
0230 end
0231 %
0232 eta_LU = curr_LU./Prf_LU;
0233 eta_LUi_1 = curr_LUi_1./Prf_LUi_1;
0234 eta_LUi_2 = curr_LUi_2./Prf_LUi_2;
0235 etaerr_LUi_1 = abs(eta_LUi_1 - eta_LU)/eta_LU;
0236 etaerr_LUi_2 = abs(eta_LUi_2 - eta_LU)/eta_LU;
0237 %
0238 disp('Incomplete LU calculation done')
0239 %
0240 % MUMPS package
0241 %
0242 if isfield(dkepath,'MUMPS') && (isfield(dkepath.MUMPS,'seq') || isfield(dkepath.MUMPS,'par')),
0243     dkeparam.invproc = -1;
0244     %
0245     [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0246     %
0247     memory_MU = dke_out.memoryLU_f_out;
0248     time_MU = dke_out.time_DKE;
0249     %
0250     nit_MU = length(dke_out.residu_f{end});
0251     %
0252     norm_MU = Znorm.x_0;
0253     curr_MU = Zcurr.x_0;
0254     Prf_MU = ZP0.x_rf_fsav;
0255     Pcoll_MU = ZP0.x_c_fsav;
0256     %
0257     disp(['Direct MUMPS calculation done'])
0258 end
0259 %
0260 % MUMPSMEX package
0261 %
0262 if exist('dmumpsmex','file'),
0263     dkeparam.invproc = -2;
0264     %
0265     [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0266     %
0267     memory_MX = dke_out.memoryLU_f_out;
0268     time_MX = dke_out.time_DKE;
0269     %
0270     nit_MX = length(dke_out.residu_f{end});
0271     %
0272     norm_MX = Znorm.x_0;
0273     curr_MX = Zcurr.x_0;
0274     Prf_MX = ZP0.x_rf_fsav;
0275     Pcoll_MX = ZP0.x_c_fsav;
0276     %
0277     disp(['MUMPSMEX calculation done'])
0278 else
0279     memory_MX = NaN;
0280     time_MX = NaN;
0281     nit_MX = NaN;
0282     norm_MX = NaN;
0283     curr_MX = NaN;
0284     Prf_MX = NaN;
0285     Pcoll_MX = NaN;
0286     %
0287     disp(['MUMPSMEX calculation skipped'])
0288 end
0289 %
0290 eta_MX = curr_MX./Prf_MX;
0291 etaerr_MX = abs(eta_MX - eta_LU)/eta_LU;
0292 %
0293 % SUPERLU package
0294 %
0295 if exist('mexlusolve'),
0296     dkeparam.invproc = -3;
0297     %
0298     [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0299     %
0300     memory_SX = dke_out.memoryLU_f_out;
0301     time_SX = dke_out.time_DKE;
0302     %
0303     nit_SX = length(dke_out.residu_f{end});
0304     %
0305     norm_SX = Znorm.x_0;
0306     curr_SX = Zcurr.x_0;
0307     Prf_SX = ZP0.x_rf_fsav;
0308     Pcoll_SX = ZP0.x_c_fsav;
0309     %
0310     disp(['SUPERLUMEX calculation done'])
0311 end
0312 %
0313 % PETSc package
0314 %
0315 if isfield(dkepath,'PETSc'),
0316     dkeparam.invproc = -30;
0317     dkeparam.PETScparam.nproc = 1;%number of processor used by PETSc (depends of the computer framework). Default is 1
0318     %
0319     [flag,message] = unix('uname -n');
0320     %
0321     if strfind(message,'saturne'),%this option runs only on the SATURNE CEA-DRFC cluster
0322         dkeparam.PETScparam.matrixtype = '-matload_type superlu_dist';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0323         dkeparam.PETScparam.kspmethod = '-ksp_type preonly';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0324         dkeparam.PETScparam.pcmethod = '-pc_type lu';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0325         %
0326         dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0327         %
0328         [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0329         %
0330         memory_PS0 = dke_out.memoryLU_f_out;
0331         time_PS0 = dke_out.time_DKE;
0332         %
0333         nit_PS0 = length(dke_out.residu_f{end});
0334         %
0335         norm_PS0 = Znorm.x_0;
0336         curr_PS0 = Zcurr.x_0;
0337         Prf_PS0 = ZP0.x_rf_fsav;
0338         Pcoll_PS0 = ZP0.x_c_fsav;
0339         %
0340         disp(['PETSc/SUPERLU calculation done (first collision operator loop in MatLab).'])
0341         %
0342         dkeparam.PETScparam.cloop = 1;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0343         %
0344         [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0345         %
0346         memory_PS1 = dke_out.memoryLU_f_out;
0347         time_PS1 = dke_out.time_DKE;
0348         %
0349         nit_PS1 = length(dke_out.residu_f{end});
0350         %
0351         norm_PS1 = Znorm.x_0;
0352         curr_PS1 = Zcurr.x_0;
0353         Prf_PS1 = ZP0.x_rf_fsav;
0354         Pcoll_PS1 = ZP0.x_c_fsav;
0355         %
0356         disp(['PETSc/SUPERLU calculation done (first collision operator loop in C).'])
0357         %
0358         flag_PETSC = 1;
0359     else
0360         flag_PETSC = 0;
0361     end
0362     %
0363     dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0364     dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0365     dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0366     %
0367     dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0368     %
0369     [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0370     %
0371     memory_PB0 = dke_out.memoryLU_f_out;
0372     time_PB0 = dke_out.time_DKE;
0373     %
0374     nit_PB0 = length(dke_out.residu_f{end});
0375     %
0376     norm_PB0 = Znorm.x_0;
0377     curr_PB0 = Zcurr.x_0;
0378     Prf_PB0 = ZP0.x_rf_fsav;
0379     Pcoll_PB0 = ZP0.x_c_fsav;
0380     %
0381     disp(['PETSc/BCGS calculation done (first collision operator loop in MatLab).'])
0382     %
0383     dkeparam.PETScparam.cloop = 1;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0384     %
0385     [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0386     %
0387     memory_PB1 = dke_out.memoryLU_f_out;
0388     time_PB1 = dke_out.time_DKE;
0389     %
0390     nit_PB1 = length(dke_out.residu_f{end});
0391     %
0392     norm_PB1 = Znorm.x_0;
0393     curr_PB1 = Zcurr.x_0;
0394     Prf_PB1 = ZP0.x_rf_fsav;
0395     Pcoll_PB1 = ZP0.x_c_fsav;
0396     %
0397     disp(['PETSc/BCGS calculation done (first collision operator loop in C).'])  
0398 end
0399 %
0400 %************************************************************************************************************************************
0401 %
0402 figure(1),clf
0403 %
0404 leg = {'CGS','BICG','QMR','MUMPSMEX'};
0405 xlim = 10.^[-9,-1];
0406 ylim = [0,2];%[0,20];%
0407 xlab = '\delta_{LU}';
0408 ylab = 'FPE calculation time (s)';
0409 tit = '';
0410 siz = 20+14i;
0411 %
0412 graph1D_jd(ludroptol_list,time_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0413 graph1D_jd(ludroptol_list,time_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0414 graph1D_jd(ludroptol_list,time_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0415 graph1D_jd(xlim,[time_MX,time_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0416 %
0417 set(gca,'ytick',[0:0.2:1]*ylim(2))
0418 set(gca,'xtick',10.^[-9:-1])
0419 set(gca,'XMinorGrid','off')
0420 set(gca,'XMinorTick','on')
0421 %
0422 % figure(2),clf
0423 % %
0424 % ylim = NaN;
0425 % ylab = 'n-1';
0426 % %
0427 % graph1D_jd(ludroptol_list,norm_LUi_1(:,1)-1,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0428 % graph1D_jd(ludroptol_list,norm_LUi_1(:,2)-1,1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0429 % graph1D_jd(ludroptol_list,norm_LUi_1(:,3)-1,1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0430 % graph1D_jd(xlim,[norm_MX-1,norm_MX-1],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0431 % %
0432 % %set(gca,'ytick',[0:0.2:1]*ylim(2))
0433 % set(gca,'xtick',10.^[-9:-1])
0434 % set(gca,'XMinorGrid','off')
0435 % set(gca,'XMinorTick','on')
0436 %
0437 figure(3),clf
0438 %
0439 ylim = [0,0.005];
0440 ylab = 'j';
0441 %
0442 graph1D_jd(ludroptol_list,curr_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0443 graph1D_jd(ludroptol_list,curr_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0444 graph1D_jd(ludroptol_list,curr_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0445 graph1D_jd(xlim,[curr_MX,curr_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0446 %
0447 set(gca,'ytick',[0:0.2:1]*ylim(2))
0448 set(gca,'xtick',10.^[-9:-1])
0449 set(gca,'XMinorGrid','off')
0450 set(gca,'XMinorTick','on')
0451 %
0452 figure(4),clf
0453 %
0454 ylim = [0,0.0002];
0455 ylab = 'P';
0456 %
0457 graph1D_jd(ludroptol_list,Prf_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0458 graph1D_jd(ludroptol_list,Prf_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0459 graph1D_jd(ludroptol_list,Prf_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0460 graph1D_jd(xlim,[Prf_MX,Prf_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0461 graph1D_jd(ludroptol_list,-Pcoll_LUi_1(:,1),1,0,'','','',NaN,xlim,ylim,'-','o','r',2,siz,gca);
0462 graph1D_jd(ludroptol_list,-Pcoll_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','o','b',2,siz,gca);
0463 graph1D_jd(ludroptol_list,-Pcoll_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','o','g',2,siz,gca);
0464 graph1D_jd(xlim,[-Pcoll_MX,-Pcoll_MX],1,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0465 %
0466 set(gca,'ytick',[0:0.2:1]*ylim(2))
0467 set(gca,'xtick',10.^[-9:-1])
0468 set(gca,'XMinorGrid','off')
0469 set(gca,'XMinorTick','on')
0470 %
0471 figure(5),clf
0472 %
0473 ylim = [0,50];
0474 ylab = '\eta';
0475 %
0476 graph1D_jd(ludroptol_list,eta_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0477 graph1D_jd(ludroptol_list,eta_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0478 graph1D_jd(ludroptol_list,eta_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0479 graph1D_jd(xlim,[eta_MX,eta_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0480 %
0481 set(gca,'ytick',[0:0.2:1]*ylim(2))
0482 set(gca,'xtick',10.^[-9:-1])
0483 set(gca,'XMinorGrid','off')
0484 set(gca,'XMinorTick','on')
0485 %
0486 figure(6),clf
0487 %
0488 ylim = 10.^[-16,0];
0489 ylab = 'Relative Error on \eta';
0490 %
0491 graph1D_jd(ludroptol_list,etaerr_LUi_1(:,1),1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0492 graph1D_jd(ludroptol_list,etaerr_LUi_1(:,2),1,1,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0493 graph1D_jd(ludroptol_list,etaerr_LUi_1(:,3),1,1,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0494 graph1D_jd(xlim,[etaerr_MX,etaerr_MX],1,1,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0495 %
0496 set(gca,'ytick',10.^[-16:2:0])
0497 set(gca,'xtick',10.^[-9:-1])
0498 set(gca,'XMinorGrid','off')
0499 set(gca,'XMinorTick','on')
0500 % %
0501 % figure(7),clf
0502 % %
0503 % ylim = [0,50];
0504 % ylab = 'Number of Legendre iterations';
0505 % %
0506 % graph1D_jd(ludroptol_list,nit_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0507 % graph1D_jd(ludroptol_list,nit_LUi_1(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0508 % graph1D_jd(ludroptol_list,nit_LUi_1(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0509 % graph1D_jd(xlim,[nit_MX,nit_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0510 % %
0511 % set(gca,'ytick',[0:0.2:1]*ylim(2))
0512 % set(gca,'xtick',10.^[-9:-1])
0513 % set(gca,'XMinorGrid','off')
0514 % set(gca,'XMinorTick','on')
0515 %
0516 figure(8),clf
0517 %
0518 leg = {'MATLAB','MUMPS'};
0519 ylim = [0,50];
0520 ylab = 'LU matrix size (MBytes)';
0521 %
0522 graph1D_jd(ludroptol_list,memory_LUi_1(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0523 graph1D_jd(xlim,[memory_MX,memory_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0524 %
0525 set(gca,'ytick',[0:0.2:1]*ylim(2))
0526 set(gca,'xtick',10.^[-9:-1])
0527 set(gca,'XMinorGrid','off')
0528 set(gca,'XMinorTick','on')
0529 %
0530 figure(9),clf
0531 %
0532 leg = {'CGS','BICG','QMR','MUMPSMEX'};
0533 xlim = 10.^[-9,-1];
0534 ylim = [0,2];%[0,20];%
0535 xlab = 'Iterative method tolerance';
0536 ylab = 'FPE calculation time (s)';
0537 tit = '';
0538 siz = 20+14i;
0539 %
0540 graph1D_jd(mlprecinv_list,time_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0541 graph1D_jd(mlprecinv_list,time_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0542 graph1D_jd(mlprecinv_list,time_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0543 graph1D_jd(xlim,[time_MX,time_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0544 %
0545 set(gca,'ytick',[0:0.2:1]*ylim(2))
0546 set(gca,'xtick',10.^[-9:-1])
0547 set(gca,'XMinorGrid','off')
0548 set(gca,'XMinorTick','on')
0549 %
0550 % figure(10),clf
0551 % %
0552 % ylim = NaN;
0553 % ylab = 'n-1';
0554 % %
0555 % graph1D_jd(mlprecinv_list,norm_LUi_2(:,1)-1,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0556 % graph1D_jd(mlprecinv_list,norm_LUi_2(:,2)-1,1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0557 % graph1D_jd(mlprecinv_list,norm_LUi_2(:,3)-1,1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0558 % graph1D_jd(xlim,[norm_MX-1,norm_MX-1],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0559 % %
0560 % %set(gca,'ytick',[0:0.2:1]*ylim(2))
0561 % set(gca,'xtick',10.^[-9:-1])
0562 % set(gca,'XMinorGrid','off')
0563 % set(gca,'XMinorTick','on')
0564 %
0565 figure(11),clf
0566 %
0567 ylim = [0,0.005];
0568 ylab = 'j';
0569 %
0570 graph1D_jd(mlprecinv_list,curr_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0571 graph1D_jd(mlprecinv_list,curr_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0572 graph1D_jd(mlprecinv_list,curr_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0573 graph1D_jd(xlim,[curr_MX,curr_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0574 %
0575 set(gca,'ytick',[0:0.2:1]*ylim(2))
0576 set(gca,'xtick',10.^[-9:-1])
0577 set(gca,'XMinorGrid','off')
0578 set(gca,'XMinorTick','on')
0579 %
0580 figure(12),clf
0581 %
0582 ylim = [0,0.0002];
0583 ylab = 'P';
0584 %
0585 graph1D_jd(mlprecinv_list,Prf_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0586 graph1D_jd(mlprecinv_list,Prf_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0587 graph1D_jd(mlprecinv_list,Prf_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0588 graph1D_jd(xlim,[Prf_MX,Prf_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0589 graph1D_jd(mlprecinv_list,-Pcoll_LUi_2(:,1),1,0,'','','',NaN,xlim,ylim,'-','o','r',2,siz,gca);
0590 graph1D_jd(mlprecinv_list,-Pcoll_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','o','b',2,siz,gca);
0591 graph1D_jd(mlprecinv_list,-Pcoll_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','o','g',2,siz,gca);
0592 graph1D_jd(xlim,[-Pcoll_MX,-Pcoll_MX],1,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0593 %
0594 set(gca,'ytick',[0:0.2:1]*ylim(2))
0595 set(gca,'xtick',10.^[-9:-1])
0596 set(gca,'XMinorGrid','off')
0597 set(gca,'XMinorTick','on')
0598 %
0599 figure(13),clf
0600 %
0601 ylim = [0,50];
0602 ylab = '\eta';
0603 %
0604 graph1D_jd(mlprecinv_list,eta_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0605 graph1D_jd(mlprecinv_list,eta_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0606 graph1D_jd(mlprecinv_list,eta_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0607 graph1D_jd(xlim,[eta_MX,eta_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0608 %
0609 set(gca,'ytick',[0:0.2:1]*ylim(2))
0610 set(gca,'xtick',10.^[-9:-1])
0611 set(gca,'XMinorGrid','off')
0612 set(gca,'XMinorTick','on')
0613 %
0614 figure(14),clf
0615 %
0616 ylim = 10.^[-12,-3];
0617 ylab = 'Relative Error on \eta';
0618 %
0619 graph1D_jd(mlprecinv_list,etaerr_LUi_2(:,1),1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0620 graph1D_jd(mlprecinv_list,etaerr_LUi_2(:,2),1,1,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0621 graph1D_jd(mlprecinv_list,etaerr_LUi_2(:,3),1,1,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0622 graph1D_jd(xlim,[etaerr_MX,etaerr_MX],1,1,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0623 %
0624 set(gca,'ytick',10.^[-12:-3])
0625 set(gca,'xtick',10.^[-9:-1])
0626 set(gca,'XMinorGrid','off')
0627 set(gca,'XMinorTick','on')
0628 %
0629 % figure(15),clf
0630 % %
0631 % ylim = [0,50];
0632 % ylab = 'Number of Legendre iterations';
0633 % %
0634 % graph1D_jd(mlprecinv_list,nit_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0635 % graph1D_jd(mlprecinv_list,nit_LUi_2(:,2),1,0,'','','',NaN,xlim,ylim,'-','+','b',2,siz,gca);
0636 % graph1D_jd(mlprecinv_list,nit_LUi_2(:,3),1,0,'','','',NaN,xlim,ylim,'-','+','g',2,siz,gca);
0637 % graph1D_jd(xlim,[nit_MX,nit_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0638 % %
0639 % set(gca,'ytick',[0:0.2:1]*ylim(2))
0640 % set(gca,'xtick',10.^[-9:-1])
0641 % set(gca,'XMinorGrid','off')
0642 % set(gca,'XMinorTick','on')
0643 %
0644 figure(16),clf
0645 %
0646 leg = {'MATLAB','MUMPS'};
0647 ylim = [0,50];
0648 ylab = 'LU matrix size (MBytes)';
0649 %
0650 graph1D_jd(mlprecinv_list,memory_LUi_2(:,1),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','+','r',2,siz,gca,0.9,0.7,0.7);
0651 graph1D_jd(xlim,[memory_MX,memory_MX],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0652 %
0653 set(gca,'ytick',[0:0.2:1]*ylim(2))
0654 set(gca,'xtick',10.^[-9:-1])
0655 set(gca,'XMinorGrid','off')
0656 set(gca,'XMinorTick','on')
0657 %
0658 print_jd(p_opt,'fig_num_lutol_time','./figures',1)
0659 %print_jd(p_opt,'fig_num_lutol_n','./figures',2)
0660 print_jd(p_opt,'fig_num_lutol_j','./figures',3)
0661 print_jd(p_opt,'fig_num_lutol_P','./figures',4)
0662 print_jd(p_opt,'fig_num_lutol_eta','./figures',5)
0663 print_jd(p_opt,'fig_num_lutol_etarel','./figures',6)
0664 %print_jd(p_opt,'fig_num_lutol_nit','./figures',7)
0665 print_jd(p_opt,'fig_num_lutol_memory','./figures',8)
0666 %
0667 print_jd(p_opt,'fig_num_ittol_time','./figures',9)
0668 %print_jd(p_opt,'fig_num_ittol_n','./figures',10)
0669 print_jd(p_opt,'fig_num_ittol_j','./figures',11)
0670 print_jd(p_opt,'fig_num_ittol_P','./figures',12)
0671 print_jd(p_opt,'fig_num_ittol_eta','./figures',13)
0672 print_jd(p_opt,'fig_num_ittol_etarel','./figures',14)
0673 %print_jd(p_opt,'fig_num_ittol_nit','./figures',15)
0674 print_jd(p_opt,'fig_num_ittol_memory','./figures',16)
0675 %
0676 delete testmatinv.txt
0677 %
0678 diary testmatinv.txt
0679 %
0680 % disp('-------------------')
0681 % disp('Direct Matlab Inversion')
0682 % disp(['droptolerance: 0 - Full CPU time: ',num2str(time_ML),' (s) - LU memory size: ',num2str(memory_ML),' (MBytes)'])
0683 % disp(['n: ',num2str(norm_ML),' - j: ',num2str(curr_ML),' - PLH: ',num2str(Prf_ML),' - Pcoll: ',num2str(Pcoll_ML)])
0684 %
0685 disp('-------------------')
0686 disp('Full LU precondtionning, CGS iterative Krylov algorithm')
0687 disp(['droptolerance: 0 - Full CPU time: ',num2str(time_LU),' (s) - LU memory size: ',num2str(memory_LU),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_LU)])
0688 disp(['n: ',num2str(norm_LU),' - j: ',num2str(curr_LU),' - PLH: ',num2str(Prf_LU),' - Pcoll: ',num2str(Pcoll_LU)])
0689 %
0690 disp('-------------------')
0691 disp('Incomplete LU precondtionning, CGS iterative Krylov algorithm')
0692 disp(['droptolerance: ',num2str(ludroptol_list(iludrop)),' - Full CPU time: ',num2str(time_LUi_1(iludrop,1)),' (s) - LU memory size: ',num2str(memory_LUi_1(iludrop,1)),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_LUi_1(iludrop,1))])
0693 disp(['n: ',num2str(norm_LUi_1(iludrop,1)),' - j: ',num2str(curr_LUi_1(iludrop,1)),' - PLH: ',num2str(Prf_LUi_1(iludrop,1)),' - Pcoll: ',num2str(Pcoll_LUi_1(iludrop,1))])
0694 %
0695 if isfield(dkepath,'MUMPS') & (isfield(dkepath.MUMPS,'seq') | isfield(dkepath.MUMPS,'par')),
0696     disp('-------------------')
0697     disp('Direct matrix inversion MUMPS')
0698     disp(['droptolerance: 0 - Full CPU time: ',num2str(time_MU),' (s) - LU memory size: ',num2str(memory_MU),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_MU)])
0699     disp(['n: ',num2str(norm_MU),' - j: ',num2str(curr_MU),' - PLH: ',num2str(Prf_MU),' - Pcoll: ',num2str(Pcoll_MU)])
0700 end
0701 if exist('dmumpsmex'),
0702     disp('-------------------')
0703     disp('Matrix inversion MUMPSMEX')
0704     disp(['droptolerance: 0 - Full CPU time: ',num2str(time_MX),' (s) - LU memory size: ',num2str(memory_MX),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_MX)])
0705     disp(['n: ',num2str(norm_MX),' - j: ',num2str(curr_MX),' - PLH: ',num2str(Prf_MX),' - Pcoll: ',num2str(Pcoll_MX)])
0706 end
0707 if exist('mexlusolve'),
0708     disp('-------------------')
0709     disp('Matrix inversion SUPERLUMEX')
0710     disp(['droptolerance: 0 - Full CPU time: ',num2str(time_SX),' (s) - LU memory size: ',num2str(memory_SX),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_SX)])
0711     disp(['n: ',num2str(norm_SX),' - j: ',num2str(curr_SX),' - PLH: ',num2str(Prf_SX),' - Pcoll: ',num2str(Pcoll_SX)])
0712 end
0713 if isfield(dkepath,'PETSc'),
0714     if flag_PETSC == 1,
0715         disp('-------------------')
0716         disp('Matrix inversion PETSc/SUPERLU, 1 processor and first order collision correction in MatLab')
0717         disp(['droptolerance: 0 - Full CPU time: ',num2str(time_PS0),' (s) - LU memory size: ',num2str(memory_PS0),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_PS0)])
0718         disp(['n: ',num2str(norm_PS0),' - j: ',num2str(curr_PS0),' - PLH: ',num2str(Prf_PS0),' - Pcoll: ',num2str(Pcoll_PS0)])
0719         disp('-------------------')
0720         disp('Matrix inversion PETSc/SUPERLU, 1 processor and first order collision correction in C')
0721         disp(['droptolerance: 0 - Full CPU time: ',num2str(time_PS1),' (s) - LU memory size: ',num2str(memory_PS1),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_PS1)])
0722         disp(['n: ',num2str(norm_PS1),' - j: ',num2str(curr_PS1),' - PLH: ',num2str(Prf_PS1),' - Pcoll: ',num2str(Pcoll_PS1)])
0723     end
0724     disp('-------------------')  
0725     disp('Matrix inversion PETSc/BCGS, 1 processor and first order collision correction in MatLab')
0726     disp(['droptolerance: 0 - Full CPU time: ',num2str(time_PB0),' (s) - LU memory size: ',num2str(memory_PB0),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_PB0)])
0727     disp(['n: ',num2str(norm_PB0),' - j: ',num2str(curr_PB0),' - PLH: ',num2str(Prf_PB0),' - Pcoll: ',num2str(Pcoll_PB0)])
0728     disp('-------------------')  
0729     disp('Matrix inversion PETSc/BCGS, 1 processor and first order collision correction in C')
0730     disp(['droptolerance: 0 - Full CPU time: ',num2str(time_PB1),' (s) - LU memory size: ',num2str(memory_PB1),' (MBytes) - Number of Leg. Iterations: ',num2str(nit_PB1)])
0731     disp(['n: ',num2str(norm_PB1),' - j: ',num2str(curr_PB1),' - PLH: ',num2str(Prf_PB1),' - Pcoll: ',num2str(Pcoll_PB1)])
0732 end
0733 %
0734 diary off
0735 %
0736 %************************************************************************************************************************************
0737 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0738 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);