run_mdce_rt

LUKE - Script for testing the distributed calculations within MatLab with call to MEX-files (using the distributed MatLab toolbox)

Script for testing the distributed calculations within MatLab with call to MEX-files (using the distributed MatLab toolbox)

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

0001 %LUKE - Script for testing the distributed calculations within MatLab with call to MEX-files (using the distributed MatLab toolbox)
0002 %
0003 %Script for testing the distributed calculations within MatLab with call to MEX-files (using the distributed MatLab toolbox)
0004 %
0005 %by Y.Peysson CEA-IRFM <yves.peysson@cea.fr> and Joan Decker CEA-IRFM (joan.decker@cea.fr)
0006 %
0007 clear all
0008 close all
0009 %
0010 dkepath = load_structures_yp('dkepath','','');
0011 %
0012 display.display_mode = 0;%no display
0013 %
0014 schedulertype = 'local';%'pbspro';%'' or 'local','xgrid',...
0015 %
0016 MDCE = 1;%Set to 1 if matlab MDCE is tested
0017 LUKEDC = 0;%Set to 1 if batch DCLUKE is tested
0018 PARLOOP = 1;%Set to 1 if parfor is tested
0019 %
0020 Nmax = Inf;%number of parallel jobs
0021 N = 12;%number of iterations
0022 %
0023 time_seq = NaN*ones(1,N);test_seq = ones(1,N);
0024 time_par = NaN*ones(1,N);test_par = NaN*ones(1,N);
0025 time_mdce = NaN*ones(1,N);test_mdce = NaN*ones(1,N);
0026 time_ldce = NaN*ones(1,N);test_ldce = NaN*ones(1,N);
0027 %
0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0029 % initial ray conditions
0030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0031 %
0032 load('EQUIL_JETliketest.mat')
0033 %
0034 omega_rf = [3.7]*2*pi*1e9;
0035 %
0036 rho0 = 0.968;
0037 theta0 = 0.0;
0038 phi0 = 0;
0039 %
0040 m0 = 0;
0041 n0 = NaN;
0042 Nphi0 = -2.0;%initial index of refraction
0043 %
0044 dNpar0 = NaN;
0045 P0_2piRp = NaN;
0046 %
0047 % Display parameters
0048 %
0049 C3POdisplay.ray = 0;
0050 C3POdisplay.equilibrium = 0;
0051 C3POdisplay.p_opt = 0;%Printing or saving option of the figures
0052 C3POdisplay.mdce = 0;%for distributed computing
0053 %
0054 % Wave parameters
0055 %
0056 waveparam.mmode = -1;%cold plasma mode [1] : (-1) m (1) p, p is the slow mode when kperp > 0 (ex : LH slow wave)
0057 waveparam.kmode = 0;%(0:cold,1:warm,2:hot;3:weak relativistic,4:full relativistic)
0058 %
0059 %Option parameter for FLR effects and cross-comparison between old FP code:
0060 %    - (0): all FLR effects
0061 %    - (1): small FLR effects and 1/vpar dependence
0062 %    - (2): small FLR effects and no 1/vpar dependence and old grid technique for DQL calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0063 %
0064 waveparam.opt_rf = NaN;
0065 %
0066 waveparam.dsmin = NaN;%minimum size for ray fragments
0067 %
0068 % -------------------------------------------------------------------------
0069 %
0070 % Global parameters for the vectorial magnetic equilibrium
0071 %
0072 fitparam.mode_equil = 1;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y)
0073 fitparam.method = 'pchip';%nearest,spline,pchip,cubic
0074 fitparam.nharm = 32;%Number of harmonics in the magnetic equilibrium interpolation (less than ntheta_equil/2)
0075 fitparam.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters
0076 fitparam.opt_load = 1;%Reload existing vectorial magnetic equilibrium (1) or overwrite it (0).
0077 %
0078 % Global parameters for the ray-tracing
0079 %
0080 rayparam.testmode = 0;
0081 rayparam.tensortype = waveparam.kmode;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0082 rayparam.t0 = 0;
0083 rayparam.tfinal = 2000;
0084 rayparam.dt0 = 1.e-4;
0085 rayparam.dS = 1.e-4;
0086 rayparam.tol = 1e-12;%when tolerance is increased (less accurate calculation of D=0), tfinal must be increased accordingly
0087 rayparam.kmax = 60000;
0088 rayparam.ncyclharm = 3;%number of cyclotron harmonics (just for hot and relativistic dielectric tensors)
0089 rayparam.reflection = 0;%1:Enforce wave reflection at plasma boundary, 0: the code calculates itself if the ray must leave of not the plasma
0090 rayparam.rel_opt = 1;%option for (1) relativistic or (0) non-relativistic calculations
0091 rayparam.nperp = 1000;%number of points in pperp integration for damping calculations
0092 rayparam.pperpmax = 10;%maximum value of pperp in damping calculations
0093 rayparam.tau_lim = 20;%value of optical depth beyond which the wave is considered absorbed
0094 %
0095 % -------------------------------------------------------------------------
0096 %
0097 % Load structures
0098 %
0099 % =========================================================================
0100 %
0101 % C3P0 ray tracing
0102 %
0103 [equil_fit] = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);%Build vectorized magnetic equilibrium structure
0104 %
0105 % Calculation of Npar0 using Nphi0
0106 %
0107 % warning : this is valid only for m=0 cases
0108 %
0109 Bx_a0_fit = ppval(equil_fit.Bx_fit.pp_a0,rho0);
0110 Bx_an_fit = ppval(equil_fit.Bx_fit.pp_an,rho0);
0111 Bx_bn_fit = ppval(equil_fit.Bx_fit.pp_bn,rho0);
0112 %
0113 By_a0_fit = ppval(equil_fit.By_fit.pp_a0,rho0);
0114 By_an_fit = ppval(equil_fit.By_fit.pp_an,rho0);
0115 By_bn_fit = ppval(equil_fit.By_fit.pp_bn,rho0);
0116 %
0117 BPHI_a0_fit = ppval(equil_fit.BPHI_fit.pp_a0,rho0);
0118 BPHI_an_fit = ppval(equil_fit.BPHI_fit.pp_an,rho0);
0119 BPHI_bn_fit = ppval(equil_fit.BPHI_fit.pp_bn,rho0);
0120 %
0121 % Build interpolated magnetic fields
0122 %
0123 [xBx] = calcval_yp(equil_fit,theta0,Bx_a0_fit,Bx_an_fit,Bx_bn_fit);
0124 [xBy] = calcval_yp(equil_fit,theta0,By_a0_fit,By_an_fit,By_bn_fit);
0125 [xBPHI] = calcval_yp(equil_fit,theta0,BPHI_a0_fit,BPHI_an_fit,BPHI_bn_fit);
0126 %
0127 xBp = sqrt(xBx.^2 + xBy.^2); 
0128 xBt = xBPHI; 
0129 xB = sqrt(xBt.^2 + xBp.^2);
0130 xBphin = xBt./xB;
0131 %
0132 Npar0 = Nphi0.*xBphin;
0133 %
0134 rayinit.omega_rf = omega_rf;
0135 %
0136 C3POparam.clustermode.main_C3PO_jd.scheduler.enforce = 0;
0137 C3POparam.clustermode.main_C3PO_jd.scheduler.display = 1;
0138 %
0139 % Loop for testing performances of the cluster
0140 %
0141 for n=1:N,
0142     %
0143     C3POparam.clustermode.main_C3PO_jd.scheduler.clustersize = min(n,Nmax);
0144     %
0145     for j = 1:n,
0146         rayinit.yNpar0 = Npar0 + 0.01*(1:n)/n;%Initial index of refraction
0147         %
0148         rayinit.yrho0 = repmat(rho0,1,n);%Initial radial position at launch
0149         rayinit.ytheta0 = repmat(theta0,1,n);%Initial poloidal position at launch
0150         rayinit.yphi0 = repmat(phi0,1,n);%Initial toroidal position at launch
0151         rayinit.ym0 = repmat(m0,1,n);%Initial poloidal mode number
0152         rayinit.yn0 = repmat(n0,1,n);%Initial toroidal mode number
0153         rayinit.ydNpar0 = repmat(dNpar0,1,n);%initial Ray spectral width
0154         rayinit.yP0_2piRp = repmat(P0_2piRp,1,n);%Lineic initial power density initial power in the ray (W/m)
0155     end
0156     %
0157     % Sequential job (for loop)
0158     %
0159     C3POparam.clustermode.main_C3PO_jd.scheduler.mode = 0;
0160     %
0161     tic;wave_seq = main_C3PO_jd(dkepath,'test',equil,equil_fit,rayinit,waveparam,fitparam,rayparam,C3POdisplay,C3POparam,[]);clear mex;clear functions;time_seq(n) = toc;
0162     %
0163     figure('Name',['Sequential calculations (for loop) for n = ',int2str(n)]),
0164     graph1D_jd(wave_seq.rays{1}.ss,wave_seq.rays{1}.srho,0,0,'s','\rho',['Sequential computing (for loop) - ',int2str(n),' cycles, dt = ',num2str(time_seq(n)),' (s)'],'',[0,max(wave_seq.rays{1}.ss)],NaN,'-','none','r',2,20,gca,0.9,0.7,0.7);
0165     for j = 2:n,
0166         graph1D_jd(wave_seq.rays{j}.ss,wave_seq.rays{j}.srho,0,0,'','','',NaN,[0,max(wave_seq.rays{1}.ss)],NaN,'-','none','r',2,20);
0167     end   
0168     %
0169     disp(['   --> time for n = ',int2str(n),' serial (for loop) calculations on (',int2str(1),' proc.) : ',num2str(time_seq(n)),' s'])
0170     disp('-------------------------------------------------------------------------')
0171     %
0172     % Matlab Distributed Computing Environment (MDCE)
0173     %
0174     if MDCE == 1,
0175         %
0176         C3POparam.clustermode.main_C3PO_jd.scheduler.mode = 1;
0177         %
0178         tic;wave_mdce = main_C3PO_jd(dkepath,'test',equil,equil_fit,rayinit,waveparam,fitparam,rayparam,C3POdisplay,C3POparam,[]);clear mex;clear functions;time_mdce(n) = toc;
0179         test_mdce(n) = comp_struct_jd(wave_seq.rays,wave_mdce.rays,0);
0180         %
0181         figure('Name',['MDCE calculations for n = ',int2str(n)]),
0182         graph1D_jd(wave_mdce.rays{1}.ss,wave_mdce.rays{1}.srho,0,0,'s','\rho',['MDCE computing - ',int2str(n),' cycles, dt = ',num2str(time_mdce(n)),' (s)'],'',[0,max(wave_mdce.rays{1}.ss)],NaN,'-','none','r',2,20,gca,0.9,0.7,0.7);
0183         for j = 2:n,
0184             graph1D_jd(wave_mdce.rays{j}.ss,wave_mdce.rays{j}.srho,0,0,'','','',NaN,[0,max(wave_mdce.rays{1}.ss)],NaN,'-','none','r',2,20);
0185         end   
0186         %
0187         disp(['   --> time for n = ',int2str(n),' MDCE calculations on (',int2str(n),' proc.) : ',num2str(time_mdce(n)),' s'])
0188         if test_mdce(n) == 1,
0189             disp('   --> sequential and MDCE calculations give same outputs');
0190         else
0191             disp('   --> WARNING: sequential and MDCE calculations give different outputs');
0192         end
0193         disp('-------------------------------------------------------------------------')
0194     end
0195     %
0196     % LUKE Distributed Computing Environment (LDCE)
0197     %
0198     if LUKEDC == 1,%All matlab may be used with this technique based on an external batch-queuing
0199         %
0200         C3POparam.clustermode.main_C3PO_jd.scheduler.mode = 2;
0201         %
0202         tic;wave_ldce = main_C3PO_jd(dkepath,'test',equil,equil_fit,rayinit,waveparam,fitparam,rayparam,C3POdisplay,C3POparam,[]);clear mex;clear functions;time_ldce(n) = toc;
0203         test_ldce(n) = comp_struct_jd(wave_seq.rays,wave_ldce.rays,0);
0204         %
0205         figure('Name',['LDCE calculations for n = ',int2str(n)]),
0206         graph1D_jd(wave_ldce.rays{1}.ss,wave_ldce.rays{1}.srho,0,0,'s','\rho',['LDCE computing - ',int2str(n),' cycles, dt = ',num2str(time_ldce(n)),' (s)'],'',[0,max(wave_ldce.rays{1}.ss)],NaN,'-','none','r',2,20,gca,0.9,0.7,0.7);
0207         for j = 2:n,
0208             graph1D_jd(wave_ldce.rays{j}.ss,wave_ldce.rays{j}.srho,0,0,'','','',NaN,[0,max(wave_ldce.rays{1}.ss)],NaN,'-','none','r',2,20);
0209         end   
0210         %
0211         disp(['   --> time for n = ',int2str(n),' LDCE calculations on (',int2str(n),' proc.) : ',num2str(time_ldce(n)),' s'])
0212         if test_ldce(n) == 1,
0213             disp('   --> sequential and LDCE calculations give same outputs');
0214         else
0215             disp('   --> WARNING: sequential and LDCE calculations give different outputs');
0216         end
0217         disp('-------------------------------------------------------------------------')
0218     end
0219     %
0220     % Parallel job (parfor loop)
0221     %
0222     if PARLOOP == 1,
0223         %
0224         C3POparam.clustermode.main_C3PO_jd.scheduler.mode = 3;
0225         %
0226         tic;wave_par = main_C3PO_jd(dkepath,'test',equil,equil_fit,rayinit,waveparam,fitparam,rayparam,C3POdisplay,C3POparam,[]);clear mex;clear functions;time_par(n) = toc;
0227         matlabpool close force
0228         test_par(n) = comp_struct_jd(wave_seq.rays,wave_par.rays,0);
0229         %
0230         figure('Name',['Parallel calculations (parfor loop) for n = ',int2str(n)]),
0231         graph1D_jd(wave_par.rays{1}.ss,wave_par.rays{1}.srho,0,0,'s','\rho',['Parallel computing (parfor loop) - ',int2str(n),' cycles, dt = ',num2str(time_par(n)),' (s)'],'',[0,max(wave_par.rays{1}.ss)],NaN,'-','none','r',2,20,gca,0.9,0.7,0.7);
0232         for j = 2:n,
0233             graph1D_jd(wave_par.rays{j}.ss,wave_par.rays{j}.srho,0,0,'','','',NaN,[0,max(wave_par.rays{1}.ss)],NaN,'-','none','r',2,20);
0234         end   
0235         %
0236         disp(['   --> time for n = ',int2str(n),' parallel (parfor loop) calculations on (',int2str(n),' proc.) : ',num2str(time_par(n)),' s'])
0237         if test_par(n) == 1,
0238             disp('   --> sequential and parallel (parfor loop) calculations give same outputs');
0239         else
0240             disp('   --> WARNING: sequential and parallel (parfor loop) calculations give different outputs');
0241         end
0242         disp('-------------------------------------------------------------------------')
0243     end
0244 end
0245 %
0246 mask_par = find(test_par);
0247 mask_mdce = find(test_mdce);
0248 mask_ldce = find(test_ldce);
0249 %
0250 figure('Name','Distributed computing'),
0251 %
0252 leg = {};
0253 %
0254 if PARLOOP == 1,
0255     leg = [leg,'Parallel'];
0256     graph1D_jd(mask_par,time_par(mask_par),0,0,'','','',NaN,[0,N+1],NaN,'none','x','r',2,20,gca);
0257 end
0258 %
0259 if MDCE == 1,
0260     leg = [leg;'MDCE'];
0261     graph1D_jd(mask_mdce,time_mdce(mask_mdce),0,0,'','','',NaN,[0,N],NaN,'none','s','g',2,20,gca);
0262 end
0263 %
0264 if LUKEDC == 1,
0265     leg = [leg;'LDCE'];
0266     graph1D_jd(mask_ldce,time_ldce(mask_ldce),0,0,'','','',NaN,[0,N],NaN,'none','d','g',2,20,gca);
0267 end
0268 %
0269 leg = ['Serial'];
0270 graph1D_jd(1:N,time_seq,0,0,'N processors','time (s)','',leg,[0,N+1],NaN,'none','o','b',2,20,gca,0.9,0.7,0.7);
0271 %