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 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_bounce';%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 = [linspace(0,0.2,11),linspace(0.25,1,16)];%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_cyl = 'TScyl';%For plasma equilibrium
0032 id_equil = 'TScirc_e1';%For plasma equilibrium
0033 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0036 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_display = 'NO_DISPLAY';%For output code display
0039 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 id_ohm = '';%For Ohmic electric contribution
0042 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0045 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_transpfaste = '';%For fast electron radial transport
0048 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 id_ripple = '';%For fast electron magnetic ripple losses
0051 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0052 %
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %************************************************************************************************************************************
0056 %
0057 [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);
0058 [equil_cyl] = load_structures_yp('equil',id_equil_cyl,path_equil);
0059 %
0060 %************************************************************************************************************************************
0061 %
0062 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0063 %Option parameter for cross-comparison between old LH code:
0064 %    - (1): 1/vpar dependence
0065 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0066 wavestruct.opt_lh = 2; % [1,1]
0067 %
0068 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0069 % or (Npar_lh,dNpar_lh) for Gaussian shape
0070 %
0071 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0072 %
0073 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0074 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0075 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0076 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0077 %
0078 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0079 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0080 %
0081 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0082 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0083 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0084 %
0085 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0086 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0087 %
0088 %************************************************************************************************************************************
0089 %
0090 if exist('dmumpsmex');dkeparam.invproc = -2;end
0091 %
0092 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)
0093 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
0094 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0095 %
0096 dkeparam.psin_S = psin_S;
0097 dkeparam.rho_S = rho_S;
0098 %
0099 % RELATIVISTIC CASE
0100 %
0101 dkeparam.pnmax_S = 30;
0102 dkeparam.np_S = 201;
0103 dkeparam.nmhu_S = 201;
0104 %
0105 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0106 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0107 %
0108 waves_cyl{1} = make_idealLHwave_jd(equil_cyl,wavestruct);
0109 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0110 %
0111 dkeparam.bounce_mode = 1;
0112 dkeparam.mhu_S = [];
0113 %
0114 [Znorm_rel,Zcurr_rel,ZP0_rel,dke_out_rel,radialDKE_rel,equilDKE_rel,momentumDKE_rel] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0115 if dke_out_rel.residu_f{end}(end) <= dkeparam.prec0_f,
0116     j_rel = Zcurr_rel.x_0_vfsav;
0117     p_rel = ZP0_rel.x_rf_fsav.';
0118     eta_rel = j_rel./p_rel;
0119 else
0120     j_rel = NaN;
0121     p_rel = NaN;
0122     eta_rel = NaN;
0123 end  
0124 %
0125 %
0126 dkeparam.bounce_mode = 0;
0127 dkeparam.mhu_S = momentumDKE_rel.mhu_S;
0128 %
0129 [Znorm_rel_cyl,Zcurr_rel_cyl,ZP0_rel_cyl,dke_out_rel_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm,waves_cyl,transpfaste,ripple,[],[]);
0130 if dke_out_rel_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0131     j_rel_cyl = Zcurr_rel_cyl.x_0;
0132     p_rel_cyl = ZP0_rel_cyl.x_rf_fsav.';
0133     eta_rel_cyl = j_rel_cyl./p_rel_cyl;
0134 else
0135     j_rel_cyl = NaN;
0136     p_rel_cyl = NaN;
0137     eta_rel_cyl = NaN;
0138 end        
0139 %
0140 % NON-RELATIVISTIC CASE
0141 %
0142 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0143 %
0144 betath = 0.001;%validated for NR limit
0145 %
0146 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0147 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0148 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0149 %
0150 equil_cyl.pTe = betath^2*mc2*ones(size(equil_cyl.pTe));
0151 equil_cyl.pzTi = betath^2*mc2*ones(size(equil_cyl.pzTi));
0152 %equil_cyl.pzTi = 1e-10*ones(size(equil_cyl.pzTi));
0153 %
0154 wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0155 wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0156 %
0157 waves_cyl{1} = make_idealLHwave_jd(equil_cyl,wavestruct);
0158 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0159 %
0160 dkeparam.pnmax_S = 10;
0161 %
0162 dkeparam.bounce_mode = 1;
0163 dkeparam.mhu_S = [];
0164 %
0165 [Znorm_nr,Zcurr_nr,ZP0_nr,dke_out_nr,radialDKE_nr,equilDKE_nr,momentumDKE_nr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0166 if dke_out_nr.residu_f{end}(end) <= dkeparam.prec0_f,
0167     j_nr = Zcurr_nr.x_0_vfsav;
0168     p_nr = ZP0_nr.x_rf_fsav.';
0169     eta_nr = j_nr./p_nr;
0170 else
0171     j_nr = NaN;
0172     p_nr = NaN;
0173     eta_nr = NaN;
0174 end  
0175 %
0176 %
0177 dkeparam.bounce_mode = 0;
0178 dkeparam.mhu_S = momentumDKE_nr.mhu_S;
0179 %
0180 [Znorm_nr_cyl,Zcurr_nr_cyl,ZP0_nr_cyl,dke_out_nr_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm,waves_cyl,transpfaste,ripple,[],[]);
0181 if dke_out_nr_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0182     j_nr_cyl = Zcurr_nr_cyl.x_0;
0183     p_nr_cyl = ZP0_nr_cyl.x_rf_fsav.';
0184     eta_nr_cyl = j_nr_cyl./p_nr_cyl;
0185 else
0186     j_nr_cyl = NaN;
0187     p_nr_cyl = NaN;
0188     eta_nr_cyl = NaN;
0189 end        
0190 %
0191 iar = equilDKE_nr.xrho*equilDKE_nr.ap/equilDKE_nr.Rp;%Circular inverse aspect ratio
0192 %
0193 %************************************************************************************************************************************
0194 %
0195 figure(1),clf
0196 %
0197 xlim = [0,1];
0198 ylim = [0,0.005];
0199 xlab = 'r/R_p';
0200 ylab = 'j';
0201 tit = '';
0202 siz = 20+14i;
0203 %
0204 graph1D_jd(iar,j_rel,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0205 graph1D_jd(iar,j_rel_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0206 %
0207 set(gca,'ytick',[0:0.2:1]*ylim(2))
0208 set(gca,'xtick',[0:0.2:1]*xlim(2))
0209 %
0210 figure(2),clf
0211 %
0212 ylim = [0,0.0005];
0213 ylab = 'P';
0214 %
0215 graph1D_jd(iar,p_rel,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0216 graph1D_jd(iar,p_rel_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0217 %
0218 set(gca,'ytick',[0:0.2:1]*ylim(2))
0219 set(gca,'xtick',[0:0.2:1]*xlim(2))
0220 %
0221 figure(3),clf
0222 %
0223 ylim = [0,50];
0224 ylab = 'j/P';
0225 %
0226 graph1D_jd(iar,eta_rel,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0227 graph1D_jd(iar,eta_rel_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0228 %
0229 set(gca,'ytick',[0:0.2:1]*ylim(2))
0230 set(gca,'xtick',[0:0.2:1]*xlim(2))
0231 %
0232 figure(4),clf
0233 %
0234 ylim = [0,0.1];
0235 ylab = 'j';
0236 %
0237 graph1D_jd(iar,j_nr,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0238 graph1D_jd(iar,j_nr_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0239 %
0240 set(gca,'ytick',[0:0.2:1]*ylim(2))
0241 set(gca,'xtick',[0:0.2:1]*xlim(2))
0242 %
0243 figure(5),clf
0244 %
0245 ylim = [0,0.01];
0246 ylab = 'P';
0247 %
0248 graph1D_jd(iar,p_nr,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0249 graph1D_jd(iar,p_nr_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0250 %
0251 set(gca,'ytick',[0:0.2:1]*ylim(2))
0252 set(gca,'xtick',[0:0.2:1]*xlim(2))
0253 %
0254 figure(6),clf
0255 %
0256 ylim = [0,20];
0257 ylab = 'j/P';
0258 %
0259 graph1D_jd(iar,eta_nr,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0260 graph1D_jd(iar,eta_nr_cyl,0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0261 %
0262 set(gca,'ytick',[0:0.2:1]*ylim(2))
0263 set(gca,'xtick',[0:0.2:1]*xlim(2))
0264 %
0265 print_jd(p_opt,'fig_LH_bounce_rel_j','./figures',1)
0266 print_jd(p_opt,'fig_LH_bounce_rel_P','./figures',2)
0267 print_jd(p_opt,'fig_LH_bounce_rel_eta','./figures',3)
0268 print_jd(p_opt,'fig_LH_bounce_nr_j','./figures',4)
0269 print_jd(p_opt,'fig_LH_bounce_nr_P','./figures',5)
0270 print_jd(p_opt,'fig_LH_bounce_nr_eta','./figures',6)
0271 %
0272 %************************************************************************************************************************************
0273 %
0274 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0275 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);