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 permission = test_permissions_yp;
0012 %
0013 if ~permission 
0014     disp('Please move the script to a local folder where you have write permission before to run it')
0015     return;
0016 end
0017 %
0018 opt_proc = 1;%test of numerical performances with number of processors at fixed Dpsipsi = 10 m2/s.
0019 nproc_max = 10;%max number of processors
0020 %
0021 mdceparam.type = 'local';%'' or 'local', 'pbspro', 'xgrid' or otherwise
0022 mdceparam.mode = 2;%distributed computing environment state: (0) disabled, (1) enabled with the dedicated Matlab toolbox which must be installed, (2) enabled with a private method (LDCE), (3) parfor parallel computing [1,1]
0023 mdceparam.enforce = 1;%enforce sequential run if distributed computing failed
0024 mdceparam.memory = 7500;%memory allocated to batch runs, in MB
0025 %
0026 p_opt = 2;
0027 %
0028 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0029 %
0030 id_simul = 'LH_karney_transp_bounce';%Simulation ID
0031 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0032 %
0033 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0034 rho_S = [0:0.025:0.65,0.7:0.1:1];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0035 %
0036 id_path = '';%For all paths used by DKE solver
0037 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0038 %
0039 id_equil = 'JETh77';%For plasma equilibrium
0040 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0041 %
0042 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0043 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0044 %
0045 id_display = 'NO_DISPLAY';%For output code display
0046 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0047 %
0048 id_ohm = '';%For Ohmic electric contribution
0049 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0050 %
0051 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0052 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0053 %
0054 id_transpfaste = 'magneticturbulence';%For fast electron radial transport
0055 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0056 %
0057 id_ripple = '';%For fast electron magnetic ripple losses
0058 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0059 %
0060 %************************************************************************************************************************************
0061 %************************************************************************************************************************************
0062 %************************************************************************************************************************************
0063 %
0064 [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);
0065 %
0066 %************************************************************************************************************************************
0067 %
0068 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0069 %Option parameter for cross-comparison between old LH code:
0070 %    - (1): 1/vpar dependence
0071 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0072 wavestruct.opt_lh = 2; % [1,1]
0073 %
0074 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0075 % or (Npar_lh,dNpar_lh) for Gaussian shape
0076 %
0077 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0078 %
0079 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0080 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0081 %
0082 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0083 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0084 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0085 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0086 %
0087 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0088 wavestruct.yD0_in_c_lh = [2];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0089 %
0090 wavestruct.yD0_in_lh_prof = [1];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0091 wavestruct.ypeak_lh = [0.4];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0092 wavestruct.ywidth_lh = [0.1];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0093 %
0094 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0095 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0096 %
0097 %************************************************************************************************************************************
0098 %
0099 if exist('dmumpsmex');dkeparam.invproc = -2;end
0100 %
0101 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)
0102 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
0103 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0104 %
0105 dkeparam.psin_S = psin_S;
0106 dkeparam.rho_S = rho_S;
0107 %
0108 %betath = 0.001;
0109 %
0110 %[qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0111 %
0112 %equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0113 %equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0114 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0115 %
0116 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0117 %
0118 % Dpsipsi scan (0 -> 100 m2/s)
0119 %
0120 Dr0_list = [0,0.001,0.01,0.1,0.3,1,3,10,30,100];
0121 Dr0_mask = [1,4,6,8];
0122 iDr0_save = 8;
0123 %
0124 nDr0 = length(Dr0_list);
0125 nr = length(rho_S) - 1;
0126 %
0127 for iDr0 = 1:nDr0,
0128     iDr0
0129     [cZcurr{iDr0},cZP0{iDr0},cdke_out{iDr0},cequilDKE{iDr0},cmksa{iDr0}] = loop_rundke(iDr0,Dr0_list,iDr0_save,id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple);
0130 end
0131 %
0132 j = NaN(nr,nDr0);
0133 P = NaN(nr,nDr0);
0134 tCPU_Dr0 = NaN(1,nDr0);
0135 Itot = NaN(1,nDr0);
0136 Ptot_2piRp = NaN(1,nDr0);
0137 %
0138 for iDr0 = 1:nDr0,
0139     %
0140     if cdke_out{iDr0}.residu_f{end}(end) <= dkeparam.prec0_f,
0141         j(:,iDr0) = cZcurr{iDr0}.x_0;
0142         %
0143         P(:,iDr0) = cZP0{iDr0}.x_rf_fsav;
0144         %
0145         tCPU_Dr0(iDr0) = cdke_out{iDr0}.time_out;
0146         %
0147         Itot(iDr0) = sum(cZcurr{iDr0}.x_0.*cequilDKE{iDr0}.xdA_dke*cmksa{iDr0}.j_ref);
0148         Ptot_2piRp(iDr0) = sum(cZP0{iDr0}.x_rf_fsav.'.*cequilDKE{iDr0}.xdV_2piRp_dke*cmksa{iDr0}.P_ref);
0149     end
0150     %
0151 end
0152 %
0153 eta = j./P;
0154 Eta = cequilDKE{1}.xne(1)*Itot./Ptot_2piRp/(2*pi)/1e19;
0155 %
0156 % N processors scan for Dpsipsi = 10 m2/s (1-> nproc_max)
0157 %
0158 if opt_proc > 0,
0159     %
0160     transpfaste.Dr0 = 10;
0161     %
0162     nproc_list = [1:nproc_max];
0163     %
0164     for inproc = 1:length(nproc_list),
0165         %
0166         [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0167         %
0168         tCPU_nproc(inproc) = dke_out.time_out;
0169         %
0170     end
0171 end
0172 %
0173 %************************************************************************************************************************************
0174 %
0175 for iDr0 = 1:length(Dr0_mask),
0176     leg{iDr0} = ['D_{\psi\psi} = ',num2str(Dr0_list(Dr0_mask(iDr0))),' m^{2}/s'];
0177 end
0178 %
0179 figure(1),clf
0180 %
0181 xlim = [0,1];
0182 ylim = [0,0.0015];
0183 xlab = 'r/a_p';
0184 ylab = 'j';
0185 tit = '';
0186 siz = 20+14i;
0187 %
0188 %for iDr0 = 1:length(Dr0_list) - 1,
0189 %    graph1D_jd(cequilDKE{1}.xrho,j(:,iDr0),i,0,0,'','','',NaN,xlim,ylim,'-','none',NaN,2,20,gca);
0190 %end
0191 graph1D_jd(cequilDKE{1}.xrho,j(:,Dr0_mask),0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',NaN,2,siz,gca,0.9,0.7,0.7);
0192 %graph1D_jd(cquilDKE{1}.xrho,besselj(0,2.4*equilDKE.xrho)*j(1,Dr0_mask(end)),0,0,'','','',NaN,xlim,ylim,'--','none','k',2,20);
0193 %
0194 set(gca,'ytick',[0:1/3:1]*ylim(2))
0195 set(gca,'xtick',[0:0.2:1]*xlim(2))
0196 %
0197 figure(2),clf
0198 %
0199 xlim = [0.01,100];
0200 ylim = [0,200];
0201 xlab = 'D_{\psi\psi} (m^{2}/s)';
0202 ylab = 'CPU time (s)';
0203 %
0204 graph1D_jd(Dr0_list(2:end),tCPU_Dr0(2:end),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0205 %
0206 set(gca,'ytick',[0:0.2:1]*ylim(2))
0207 set(gca,'xtick',[0.01,0.1,1,10,100])
0208 set(gca,'XMinorGrid','off')
0209 set(gca,'XMinorTick','on')
0210 %
0211 figure(3),clf
0212 %
0213 ylim = [0,0.06];
0214 ylab = 'I (MA)';
0215 %
0216 graph1D_jd(Dr0_list(2:end),-Itot(2:end),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0217 %
0218 set(gca,'ytick',[0:0.25:1]*ylim(2))
0219 set(gca,'xtick',[0.01,0.1,1,10,100])
0220 set(gca,'XMinorGrid','off')
0221 set(gca,'XMinorTick','on')
0222 %
0223 figure(4),clf
0224 %
0225 ylim = [0,0.004];
0226 ylab = 'P/(2\piR_p) (MW/m)';
0227 %
0228 graph1D_jd(Dr0_list(2:end),Ptot_2piRp(2:end),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0229 %
0230 set(gca,'ytick',[0:0.25:1]*ylim(2))
0231 set(gca,'xtick',[0.01,0.1,1,10,100])
0232 set(gca,'XMinorGrid','off')
0233 set(gca,'XMinorTick','on')
0234 %
0235 figure(5),clf
0236 %
0237 ylim = [0,6];
0238 ylab = 'n_{19}IR_p/P';
0239 %
0240 graph1D_jd(Dr0_list(2:end),-Eta(2:end),1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0241 %
0242 set(gca,'ytick',[0:0.25:1]*ylim(2))
0243 set(gca,'xtick',[0.01,0.1,1,10,100])
0244 set(gca,'XMinorGrid','off')
0245 set(gca,'XMinorTick','on')
0246 %
0247 print_jd(p_opt,'fig_JET_transp_nobounce_j')
0248 print_jd(p_opt,'fig_JET_transp_nobounce_cpu')
0249 print_jd(p_opt,'fig_JET_transp_nobounce_I')
0250 print_jd(p_opt,'fig_JET_transp_nobounce_P')
0251 print_jd(p_opt,'fig_JET_transp_nobounce_Eta')
0252 %
0253 if opt_proc > 0, 
0254     %
0255     figure(6),clf
0256     %
0257     xlim = [0,length(nproc_list)+1];
0258     ylim = [0,round(max(tCPU_nproc)*2)];
0259     xlab = 'Number of processors';
0260     ylab = 'CPU time (s)';
0261     tit = 'D_{\psi\psi} = 10 (m^{2}/s)';
0262     %
0263     graph1D_jd(nproc_list,tCPU_nproc,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0264     %
0265     print_jd(p_opt,'fig_transp_nobounce_cpu_proc')
0266     %
0267 end
0268 %
0269 %************************************************************************************************************************************
0270 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0271 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);