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_Zi';%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 = 201; 0098 dkeparam.nmhu_S = 201; 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 Zi_list = [1:100];%Z = 74 is tungsten fully stripped ion 0112 % 0113 nZi = length(Zi_list); 0114 % 0115 j_0 = NaN(1,nZi); 0116 j_1 = NaN(1,nZi); 0117 j_2 = NaN(1,nZi); 0118 %j_4 = NaN(1,nZi); 0119 % 0120 P_0 = NaN(1,nZi); 0121 P_1 = NaN(1,nZi); 0122 P_2 = NaN(1,nZi); 0123 %P_4 = NaN(1,nZi); 0124 % 0125 for iZi = 1:nZi, 0126 % 0127 Zi = Zi_list(iZi); 0128 equil.zZi = [1,1,1,Zi]; 0129 equil.pzni = [0;0;0;1/Zi]*equil.pne; 0130 % 0131 waves{1} = make_idealLHwave_jd(equil,wavestruct); 0132 % 0133 dkeparam.coll_mode = 0;% Relativistic Maxwellian background 0134 [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0135 j_0(iZi) = Zcurr.x_0; 0136 P_0(iZi) = ZP0.x_rf_fsav; 0137 % 0138 dkeparam.coll_mode = 1;% High-velocity limit 0139 [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0140 j_1(iZi) = Zcurr.x_0; 0141 P_1(iZi) = ZP0.x_rf_fsav; 0142 % 0143 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker 0144 [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0145 if dke_out.residu_f{end}(end) <= dkeparam.prec0_f, 0146 j_2(iZi) = Zcurr.x_0; 0147 P_2(iZi) = ZP0.x_rf_fsav; 0148 end 0149 % 0150 % dkeparam.coll_mode = 4;% Linearized Belaiev-Budker 0151 % [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0152 % j_4(iZi) = Zcurr.x_0; 0153 % P_4(iZi) = ZP0.x_rf_fsav; 0154 % 0155 end 0156 % 0157 eta_0 = j_0./P_0; 0158 eta_1 = j_1./P_1; 0159 eta_2 = j_2./P_2; 0160 %eta_4 = j_4./P_4; 0161 % 0162 j_0_nr_Karney_1 = 0.05759; 0163 P_0_nr_Karney_1 = 0.004012; 0164 eta_0_nr_Karney_1 = 14.35; 0165 % 0166 j_2_nr_Karney_1 = 0.07092; 0167 P_2_nr_Karney_1 = 0.004294; 0168 eta_2_nr_Karney_1 = 16.52; 0169 % 0170 %eta_0_nr_Karney = eta_0_nr_Karney_1*6./(5 + Zi_list); 0171 eta_2_nr_Karney = eta_2_nr_Karney_1*6./(5 + Zi_list); 0172 % 0173 %************************************************************************************************************************************ 0174 % 0175 figure(1),clf 0176 % 0177 %leg = {'Linearized','High v limit','Maxwellian NR','Maxwellian'}; 0178 leg = {'Linearized','High v limit','Maxwellian'}; 0179 xlim = [0,max(Zi_list)]; 0180 ylim = [0,0.1]; 0181 xlab = 'Z_i'; 0182 ylab = 'j'; 0183 tit = ''; 0184 siz = 20+14i; 0185 % 0186 graph1D_jd(Zi_list,j_2,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7); 0187 graph1D_jd(Zi_list,j_1,0,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca); 0188 %graph1D_jd(Zi_list,j_4,0,0,'','','',NaN,xlim,ylim,'-','none','m',2,siz,gca); 0189 graph1D_jd(Zi_list,j_0,0,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca); 0190 %graph1D_jd(xlim,[j_2_nr_Karney_1,j_2_nr_Karney_1],0,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca); 0191 %graph1D_jd(xlim,[j_0_nr_Karney_1,j_0_nr_Karney_1],0,0,'','','',NaN,xlim,ylim,'--','none','m',2,siz,gca); 0192 % 0193 set(gca,'ytick',[0:0.2:1]*ylim(2)) 0194 set(gca,'xtick',[0:0.2:1]*xlim(2)) 0195 % 0196 figure(2),clf 0197 % 0198 ylim = [0,0.02]; 0199 ylab = 'P'; 0200 % 0201 graph1D_jd(Zi_list,P_2,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7); 0202 graph1D_jd(Zi_list,P_1,0,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca); 0203 %graph1D_jd(Zi_list,P_4,0,0,'','','',NaN,xlim,ylim,'-','none','m',2,siz,gca); 0204 graph1D_jd(Zi_list,P_0,0,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca); 0205 %graph1D_jd(xlim,[P_2_nr_Karney_1,P_2_nr_Karney_1],0,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca); 0206 %graph1D_jd(xlim,[P_0_nr_Karney_1,P_0_nr_Karney_1],0,0,'','','',NaN,xlim,ylim,'--','none','m',2,siz,gca); 0207 % 0208 set(gca,'ytick',[0:0.2:1]*ylim(2)) 0209 set(gca,'xtick',[0:0.2:1]*xlim(2)) 0210 % 0211 figure(3),clf 0212 % 0213 ylim = [0,20]; 0214 ylab = 'j/P'; 0215 % 0216 graph1D_jd(Zi_list,eta_2,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7); 0217 graph1D_jd(Zi_list,eta_1,0,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca); 0218 %graph1D_jd(Zi_list,eta_4,0,0,'','','',NaN,xlim,ylim,'-','none','m',2,siz,gca); 0219 graph1D_jd(Zi_list,eta_0,0,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca); 0220 graph1D_jd(Zi_list,eta_2_nr_Karney,0,0,'','','',NaN,xlim,ylim,'--','none','k',0.5,siz,gca); 0221 %graph1D_jd(Zi_list,eta_0_nr_Karney,0,0,'','','',NaN,xlim,ylim,'--','none','m',2,siz,gca); 0222 % 0223 set(gca,'ytick',[0:0.2:1]*ylim(2)) 0224 set(gca,'xtick',[0:0.2:1]*xlim(2)) 0225 % 0226 print_jd(p_opt,'fig_j_Zi','./figures',1) 0227 print_jd(p_opt,'fig_P_Zi','./figures',2) 0228 print_jd(p_opt,'fig_eta_Zi','./figures',3) 0229 % 0230 %************************************************************************************************************************************ 0231 % 0232 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']); 0233 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);