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