rundke_D0

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_D0';%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 = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0072 wavestruct.yvparmax_lh = [5];%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 if exist('dmumpsmex');dkeparam.invproc = -2;end
0092 %
0093 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)
0094 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
0095 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0096 %
0097 dkeparam.np_S = 401;
0098 dkeparam.nmhu_S = 401;
0099 dkeparam.pnmax_S = 10;
0100 %
0101 dkeparam.psin_S = psin_S;
0102 dkeparam.rho_S = rho_S;
0103 %
0104 betath = 0.001;
0105 %
0106 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0107 %
0108 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0109 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0110 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0111 %
0112 D0_list = [logspace(-6,-2,9),logspace(-1.75,2,16)];
0113 %
0114 nD0 = length(D0_list);
0115 %
0116 j_0 = NaN(1,nD0);
0117 j_1 = NaN(1,nD0);
0118 j_2 = NaN(1,nD0);
0119 %
0120 j_0a = NaN(1,nD0);
0121 j_1a = NaN(1,nD0);
0122 j_2a = NaN(1,nD0);
0123 %
0124 P_0 = NaN(1,nD0);
0125 P_1 = NaN(1,nD0);
0126 P_2 = NaN(1,nD0);
0127 %
0128 P_0a = NaN(1,nD0);
0129 P_1a = NaN(1,nD0);
0130 P_2a = NaN(1,nD0);
0131 %
0132 Pc_0 = NaN(1,nD0);
0133 Pc_1 = NaN(1,nD0);
0134 Pc_2 = NaN(1,nD0);
0135 %
0136 Pc_0a = NaN(1,nD0);
0137 Pc_1a = NaN(1,nD0);
0138 Pc_2a = NaN(1,nD0);
0139 %
0140 for iD0 = 1:nD0,
0141     %
0142     wavestruct.yD0_in_c_lh = [D0_list(iD0)];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0143     %
0144     waves{1} = make_idealLHwave_jd(equil,wavestruct);
0145     %
0146     dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0147     %
0148     dkeparam.oldcross_mode = 0;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0149     [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0150     j_0(iD0) = Zcurr.x_0;
0151     P_0(iD0) = ZP0.x_rf_fsav;
0152     Pc_0(iD0) = -ZP0.x_c_fsav;
0153     %
0154     dkeparam.oldcross_mode = 1;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0155     [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0156     j_0a(iD0) = Zcurr.x_0;
0157     P_0a(iD0) = ZP0.x_rf_fsav;
0158     Pc_0a(iD0) = -ZP0.x_c_fsav;
0159     %
0160     dkeparam.coll_mode = 1;% High-velocity limit
0161     %
0162     dkeparam.oldcross_mode = 0;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0163     [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0164     j_1(iD0) = Zcurr.x_0;
0165     P_1(iD0) = ZP0.x_rf_fsav;
0166     Pc_1(iD0) = -ZP0.x_c_fsav;
0167     %
0168     dkeparam.oldcross_mode = 1;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0169     [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0170     j_1a(iD0) = Zcurr.x_0;
0171     P_1a(iD0) = ZP0.x_rf_fsav;
0172     Pc_1a(iD0) = -ZP0.x_c_fsav;
0173     %
0174     dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0175     %
0176     dkeparam.oldcross_mode = 0;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0177     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0178     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0179         j_2(iD0) = Zcurr.x_0;
0180         P_2(iD0) = ZP0.x_rf_fsav;
0181         Pc_2(iD0) = -ZP0.x_c_fsav;
0182     end
0183     %
0184     dkeparam.oldcross_mode = 1;%Cross-derivatives as given in paper Peysson-Shoucri CPC 98
0185     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0186     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0187         j_2a(iD0) = Zcurr.x_0;
0188         P_2a(iD0) = ZP0.x_rf_fsav;
0189         Pc_2a(iD0) = -ZP0.x_c_fsav; 
0190     end
0191     %
0192 end
0193 %
0194 eta_0 = j_0./P_0;
0195 eta_1 = j_1./P_1;
0196 eta_2 = j_2./P_2;
0197 %
0198 eta_0a = j_0a./P_0a;
0199 eta_1a = j_1a./P_1a;
0200 eta_2a = j_2a./P_2a;
0201 %
0202 %************************************************************************************************************************************
0203 %
0204 figure(1),clf
0205 %
0206 %leg = {'Linearized','High v limit','Maxwellian NR','Maxwellian'};
0207 leg = {'Linearized','High v limit','Maxwellian'};
0208 xlim = 10.^[-5,1];
0209 xticks = [1e-06 1e-5 0.0001 0.001 0.01 0.1 1 10 100];
0210 ylim = 10.^[-6,-1];
0211 xlab = 'D_0';
0212 ylab = 'j';
0213 tit = '';
0214 siz = 20+14i;
0215 %
0216 graph1D_jd(D0_list,j_2,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0217 graph1D_jd(D0_list,j_1,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0218 graph1D_jd(D0_list,j_0,1,1,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0219 %
0220 set(gca,'ytick',[1e-06 1e-05 0.0001 0.001 0.01 0.1])
0221 set(gca,'xtick',xticks)
0222 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0223 set(gca,'XMinorGrid','off')
0224 set(gca,'XMinorTick','on')
0225 set(gca,'YMinorGrid','off')
0226 set(gca,'YMinorTick','on')
0227 %
0228 figure(2),clf
0229 %
0230 ylim = 10.^[-7,-2];
0231 ylab = 'P';
0232 %
0233 graph1D_jd(D0_list,P_2,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0234 graph1D_jd(D0_list,P_1,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0235 graph1D_jd(D0_list,P_0,1,1,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0236 graph1D_jd(D0_list,Pc_2,1,1,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0237 graph1D_jd(D0_list,Pc_1,1,1,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0238 graph1D_jd(D0_list,Pc_0,1,1,'','','',leg,xlim,ylim,'--','none','g',2,siz,gca);
0239 %
0240 set(gca,'ytick',[1e-07 1e-06 1e-05 0.0001 0.001 0.01])
0241 set(gca,'xtick',xticks)
0242 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0243 set(gca,'XMinorGrid','off')
0244 set(gca,'XMinorTick','on')
0245 set(gca,'YMinorGrid','off')
0246 set(gca,'YMinorTick','on')
0247 %
0248 figure(3),clf
0249 %
0250 ylim = [0,20];
0251 ylab = 'j/P';
0252 %
0253 graph1D_jd(D0_list,eta_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0254 graph1D_jd(D0_list,eta_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0255 graph1D_jd(D0_list,eta_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0256 %
0257 set(gca,'ytick',[0:0.2:1]*ylim(2))
0258 set(gca,'xtick',xticks)
0259 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0260 set(gca,'XMinorGrid','off')
0261 set(gca,'XMinorTick','on')
0262 %
0263 figure(4),clf
0264 %
0265 %leg = {'Linearized','High v limit','Maxwellian NR','Maxwellian'};
0266 leg = {'Linearized','High v limit','Maxwellian'};
0267 xlim = 10.^[-5,1];
0268 xticks = [1e-06 1e-5 0.0001 0.001 0.01 0.1 1 10 100];
0269 ylim = 10.^[-6,-1];
0270 xlab = 'D_0';
0271 ylab = 'j';
0272 tit = '';
0273 %
0274 graph1D_jd(D0_list,j_2a,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0275 graph1D_jd(D0_list,j_1a,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0276 graph1D_jd(D0_list,j_0a,1,1,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0277 %
0278 set(gca,'ytick',[1e-06 1e-05 0.0001 0.001 0.01 0.1])
0279 set(gca,'xtick',xticks)
0280 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0281 set(gca,'XMinorGrid','off')
0282 set(gca,'XMinorTick','on')
0283 set(gca,'YMinorGrid','off')
0284 set(gca,'YMinorTick','on')
0285 %
0286 figure(5),clf
0287 %
0288 ylim = 10.^[-7,-2];
0289 ylab = 'P';
0290 %
0291 graph1D_jd(D0_list,P_2a,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0292 graph1D_jd(D0_list,P_1a,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0293 graph1D_jd(D0_list,P_0a,1,1,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0294 graph1D_jd(D0_list,Pc_2a,1,1,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0295 graph1D_jd(D0_list,Pc_1a,1,1,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0296 graph1D_jd(D0_list,Pc_0a,1,1,'','','',leg,xlim,ylim,'--','none','g',2,siz,gca);
0297 %
0298 set(gca,'ytick',[1e-07 1e-06 1e-05 0.0001 0.001 0.01])
0299 set(gca,'xtick',xticks)
0300 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0301 set(gca,'XMinorGrid','off')
0302 set(gca,'XMinorTick','on')
0303 set(gca,'YMinorGrid','off')
0304 set(gca,'YMinorTick','on')
0305 %
0306 figure(6),clf
0307 %
0308 ylim = [0,20];
0309 ylab = 'j/P';
0310 %
0311 graph1D_jd(D0_list,eta_2a,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0312 graph1D_jd(D0_list,eta_1a,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0313 graph1D_jd(D0_list,eta_0a,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0314 %
0315 set(gca,'ytick',[0:0.2:1]*ylim(2))
0316 set(gca,'xtick',xticks)
0317 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0318 set(gca,'XMinorGrid','off')
0319 set(gca,'XMinorTick','on')
0320 %
0321 print_jd(p_opt,'fig_j_D0','./figures',1)
0322 print_jd(p_opt,'fig_P_D0','./figures',2)
0323 print_jd(p_opt,'fig_eta_D0','./figures',3)
0324 %
0325 print_jd(p_opt,'fig_j_D0_oldcross','./figures',4)
0326 print_jd(p_opt,'fig_P_D0_oldcross','./figures',5)
0327 print_jd(p_opt,'fig_eta_D0_oldcross','./figures',6)
0328 %
0329 %************************************************************************************************************************************
0330 %
0331 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0332 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);