rundke_time

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 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0014 %
0015 id_simul = 'LH_karney_time';%Simulation ID
0016 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0017 %
0018 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0019 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)
0020 %
0021 id_path = '';%For all paths used by DKE solver
0022 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0023 %
0024 id_equil = 'TScyl';%For plasma equilibrium
0025 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0026 %
0027 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0028 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0029 %
0030 id_display = 'NO_DISPLAY';%For output code display
0031 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0032 %
0033 id_ohm = '';%For Ohmic electric contribution
0034 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0035 %
0036 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0037 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0038 %
0039 id_transpfaste = '';%For fast electron radial transport
0040 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0041 %
0042 id_ripple = '';%For fast electron magnetic ripple losses
0043 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0044 %
0045 %************************************************************************************************************************************
0046 %************************************************************************************************************************************
0047 %************************************************************************************************************************************
0048 %
0049 [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);
0050 %
0051 %************************************************************************************************************************************
0052 %
0053 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0054 %Option parameter for cross-comparison between old LH code:
0055 %    - (1): 1/vpar dependence
0056 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0057 wavestruct.opt_lh = 2; % [1,1]
0058 %
0059 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0060 % or (Npar_lh,dNpar_lh) for Gaussian shape
0061 %
0062 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0063 %
0064 wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0065 wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0066 %
0067 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0068 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0069 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0070 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0071 %
0072 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0073 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0074 %
0075 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0076 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0077 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0078 %
0079 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0080 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0081 %
0082 %************************************************************************************************************************************
0083 %
0084 if exist('dmumpsmex');dkeparam.invproc = -2;end
0085 %
0086 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)
0087 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
0088 dkeparam.timevol = 1;%to calculate moments at all internal times
0089 %
0090 dkeparam.nmhu_S = 201;
0091 dkeparam.np_S = 201;
0092 %
0093 dkeparam.psin_S = psin_S;
0094 dkeparam.rho_S = rho_S;
0095 %
0096 nit = 141;
0097 tn = logspace(-2,5,nit);
0098 %
0099 dkeparam.tn = tn;
0100 dkeparam.dtn = NaN;%specified by tn vector
0101 %
0102 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0103 %
0104 betath = 0.001;%validated for NR limit
0105 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0106 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0107 %
0108 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0109 %
0110 dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0111 [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0112 %
0113 j_0 = NaN(1,nit);
0114 P_0 = NaN(1,nit);
0115 for it = 1:nit,
0116     j_0(it) = Zcurr(it).x_0;
0117     P_0(it) = ZP0(it).x_rf_fsav;
0118 end
0119 %
0120 dkeparam.coll_mode = 1;% High-velocity limit
0121 [dummy,Zcurr,ZP0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0122 %
0123 j_1 = NaN(1,nit);
0124 P_1 = NaN(1,nit);
0125 for it = 1:nit,
0126     j_1(it) = Zcurr(it).x_0;
0127     P_1(it) = ZP0(it).x_rf_fsav;
0128 end
0129 %
0130 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0131 [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0132 %
0133 residue_2 = NaN(1,nit);
0134 j_2 = NaN(1,nit);
0135 P_2 = NaN(1,nit);
0136 for it = 1:nit,
0137     residue_2(it) = dke_out.residu_f{it}(end);
0138     j_2(it) = Zcurr(it).x_0;
0139     P_2(it) = ZP0(it).x_rf_fsav;
0140 end
0141 j_2(residue_2 > dkeparam.prec0_f) = NaN;
0142 P_2(residue_2 > dkeparam.prec0_f) = NaN;
0143 %
0144 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)
0145 dkeparam.norm_mode_f = 0;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0146 %
0147 [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0148 %
0149 residue_2_0 = NaN(1,nit);
0150 j_2_0 = NaN(1,nit);
0151 P_2_0 = NaN(1,nit);
0152 for it = 1:nit,
0153     residue_2_0(it) = dke_out.residu_f{it}(end);
0154     j_2_0(it) = Zcurr(it).x_0;
0155     P_2_0(it) = ZP0(it).x_rf_fsav;
0156 end
0157 j_2_0(residue_2_0 > dkeparam.prec0_f) = NaN;
0158 P_2_0(residue_2_0 > dkeparam.prec0_f) = NaN;
0159 %
0160 dkeparam.boundary_mode_f = 1;%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)
0161 dkeparam.norm_mode_f = 0;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0162 %
0163 [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0164 %
0165 residue_2_1 = NaN(1,nit);
0166 j_2_1 = NaN(1,nit);
0167 P_2_1 = NaN(1,nit);
0168 for it = 1:nit,
0169     residue_2_1(it) = dke_out.residu_f{it}(end);
0170     j_2_1(it) = Zcurr(it).x_0;
0171     P_2_1(it) = ZP0(it).x_rf_fsav;
0172 end
0173 j_2_1(residue_2_1 > dkeparam.prec0_f) = NaN;
0174 P_2_1(residue_2_1 > dkeparam.prec0_f) = NaN;
0175 %
0176 eta_0 = j_0./P_0;
0177 eta_1 = j_1./P_1;
0178 eta_2 = j_2./P_2;
0179 eta_2_0 = j_2_0./P_2_0;
0180 eta_2_1 = j_2_1./P_2_1;
0181 %
0182 j_0_nr_Karney = 0.05759;
0183 P_0_nr_Karney = 0.004012;
0184 eta_0_nr_Karney = 14.35;
0185 %
0186 j_2_nr_Karney = 0.07092;
0187 P_2_nr_Karney = 0.004294;
0188 eta_2_nr_Karney = 16.52;
0189 %
0190 %************************************************************************************************************************************
0191 %
0192 %
0193 figure(1),clf
0194 %
0195 leg = {'Linearized','High v limit','Maxwellian'};
0196 xlim = 10.^[-2,5];
0197 ylim = [0,0.1];
0198 xlab = 't/\tau_c';
0199 ylab = 'j';
0200 tit = '';
0201 siz = 20+14i;
0202 %
0203 graph1D_jd(tn,j_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0204 graph1D_jd(tn,j_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0205 graph1D_jd(tn,j_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0206 graph1D_jd(xlim,[j_2_nr_Karney,j_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0207 graph1D_jd(xlim,[j_0_nr_Karney,j_0_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','g',2,siz,gca);
0208 %
0209 set(gca,'ytick',[0:0.02:0.1])
0210 set(gca,'xtick',10.^[-2:5])
0211 set(gca,'XMinorGrid','off')
0212 set(gca,'XMinorTick','on')
0213 %
0214 figure(2),clf
0215 %
0216 ylim = [0,0.015];
0217 ylab = 'P';
0218 %
0219 graph1D_jd(tn,P_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0220 graph1D_jd(tn,P_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0221 graph1D_jd(tn,P_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0222 graph1D_jd(xlim,[P_2_nr_Karney,P_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0223 graph1D_jd(xlim,[P_0_nr_Karney,P_0_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','g',2,siz,gca);
0224 %
0225 set(gca,'ytick',[0:0.0025:0.015])
0226 set(gca,'xtick',10.^[-2:5])
0227 set(gca,'XMinorGrid','off')
0228 set(gca,'XMinorTick','on')
0229 %
0230 figure(3),clf
0231 %
0232 ylim = [0,20];
0233 ylab = 'j/P';
0234 %
0235 graph1D_jd(tn,eta_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0236 graph1D_jd(tn,eta_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0237 graph1D_jd(tn,eta_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0238 graph1D_jd(xlim,[eta_2_nr_Karney,eta_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0239 graph1D_jd(xlim,[eta_0_nr_Karney,eta_0_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','g',2,siz,gca);
0240 %
0241 set(gca,'ytick',[0:4:20])
0242 set(gca,'xtick',10.^[-2:5])
0243 set(gca,'XMinorGrid','off')
0244 set(gca,'XMinorTick','on')
0245 %
0246 figure(4),clf
0247 %
0248 leg = {'Free Mode','Fixed point','Normalization'};
0249 ylim = [0,0.1];
0250 ylab = 'j';
0251 %
0252 graph1D_jd(tn,j_2_0,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','b',2,siz,gca,0.9,0.7,0.7);
0253 graph1D_jd(tn,j_2_1,1,0,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0254 graph1D_jd(tn,j_2,1,0,'','','',leg,xlim,ylim,'-','none','r',2,siz,gca);
0255 graph1D_jd(xlim,[j_2_nr_Karney,j_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0256 %
0257 set(gca,'ytick',[0:0.02:0.1])
0258 set(gca,'xtick',10.^[-2:5])
0259 set(gca,'XMinorGrid','off')
0260 set(gca,'XMinorTick','on')
0261 %
0262 figure(5),clf
0263 %
0264 ylim = [0,0.015];
0265 ylab = 'P';
0266 %
0267 graph1D_jd(tn,P_2_0,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','b',2,siz,gca,0.9,0.7,0.7);
0268 graph1D_jd(tn,P_2_1,1,0,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0269 graph1D_jd(tn,P_2,1,0,'','','',leg,xlim,ylim,'-','none','r',2,siz,gca);
0270 graph1D_jd(xlim,[P_2_nr_Karney,P_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0271 %
0272 set(gca,'ytick',[0:0.0025:0.015])
0273 set(gca,'xtick',10.^[-2:5])
0274 set(gca,'XMinorGrid','off')
0275 set(gca,'XMinorTick','on')
0276 %
0277 figure(6),clf
0278 %
0279 ylim = [0,20];
0280 ylab = 'j/P';
0281 %
0282 graph1D_jd(tn,eta_2_0,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','b',2,siz,gca,0.9,0.7,0.7);
0283 graph1D_jd(tn,eta_2_1,1,0,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0284 graph1D_jd(tn,eta_2,1,0,'','','',leg,xlim,ylim,'-','none','r',2,siz,gca);
0285 graph1D_jd(xlim,[eta_2_nr_Karney,eta_2_nr_Karney],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0286 %
0287 set(gca,'ytick',[0:4:20])
0288 set(gca,'xtick',10.^[-2:5])
0289 set(gca,'XMinorGrid','off')
0290 set(gca,'XMinorTick','on')
0291 %
0292 print_jd(p_opt,'fig_j_coll_time','./figures',1)
0293 print_jd(p_opt,'fig_P_coll_time','./figures',2)
0294 print_jd(p_opt,'fig_eta_coll_time','./figures',3)
0295 print_jd(p_opt,'fig_j_norm_time','./figures',4)
0296 print_jd(p_opt,'fig_P_norm_time','./figures',5)
0297 print_jd(p_opt,'fig_eta_norm_time','./figures',6)
0298 %
0299 %************************************************************************************************************************************
0300 %
0301 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0302 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);