rundke_grid

PURPOSE ^

Script for running the DKE solver (can be modified by the user for specific simulations)

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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_grid';%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.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0072 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0073 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0074 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0075 %
0076 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0077 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0078 %
0079 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0080 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0081 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0082 %
0083 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0084 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0085 %
0086 %************************************************************************************************************************************
0087 %
0088 if exist('dmumpsmex');dkeparam.invproc = -2;end
0089 %
0090 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)
0091 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
0092 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0093 %
0094 dkeparam.psin_S = psin_S;
0095 dkeparam.rho_S = rho_S;
0096 %
0097 np_list = round(logspace(1,3,41));
0098 %
0099 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0100 %
0101 nnp = length(np_list);
0102 %
0103 j_rel_50 = NaN(1,nnp);
0104 j_rel_100 = NaN(1,nnp);
0105 j_rel_200 = NaN(1,nnp);
0106 p_rel_50 = NaN(1,nnp);
0107 p_rel_100 = NaN(1,nnp);
0108 p_rel_200 = NaN(1,nnp);
0109 %
0110 j_nr_50 = NaN(1,nnp);
0111 j_nr_100 = NaN(1,nnp);
0112 j_nr_200 = NaN(1,nnp);
0113 p_nr_50 = NaN(1,nnp);
0114 p_nr_100 = NaN(1,nnp);
0115 p_nr_200 = NaN(1,nnp);
0116 %
0117 for inp = 1:nnp,
0118     %
0119     dkeparam.np_S = np_list(inp) + 1;
0120     %
0121     % RELATIVISTIC CASE
0122     %
0123     betath = 0.1;
0124     equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0125     equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0126     %
0127     wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0128     wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0129     %
0130     waves{1} = make_idealLHwave_jd(equil,wavestruct);
0131     %
0132     dkeparam.pnmax_S = 30;
0133     %
0134     dkeparam.nmhu_S = 51;
0135     %
0136     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0137     %
0138     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0139         j_rel_50(inp) = Zcurr.x_0;
0140         p_rel_50(inp) = ZP0.x_rf_fsav;
0141     end
0142     %
0143     dkeparam.nmhu_S = 101;
0144     %
0145     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0146     %
0147     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0148         j_rel_100(inp) = Zcurr.x_0;
0149         p_rel_100(inp) = ZP0.x_rf_fsav;
0150     end
0151     %
0152     dkeparam.nmhu_S = 201;
0153     %
0154     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0155     %
0156     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0157         j_rel_200(inp) = Zcurr.x_0;
0158         p_rel_200(inp) = ZP0.x_rf_fsav;
0159     end
0160     %
0161     % NON-RELATIVISTIC CASE
0162     %
0163     betath = 0.001;
0164     equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0165     equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0166     %
0167     wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0168     wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0169     %
0170     waves{1} = make_idealLHwave_jd(equil,wavestruct);
0171     %
0172     dkeparam.pnmax_S = 20;
0173     %
0174     dkeparam.nmhu_S = 51;
0175     %
0176     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0177     %
0178     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0179         j_nr_50(inp) = Zcurr.x_0;
0180         p_nr_50(inp) = ZP0.x_rf_fsav;
0181     end
0182     %
0183     dkeparam.nmhu_S = 101;
0184     %
0185     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0186     %
0187     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0188         j_nr_100(inp) = Zcurr.x_0;
0189         p_nr_100(inp) = ZP0.x_rf_fsav;
0190     end
0191     %
0192     dkeparam.nmhu_S = 201;
0193     %
0194     [dummy,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0195     %
0196     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0197         j_nr_200(inp) = Zcurr.x_0;
0198         p_nr_200(inp) = ZP0.x_rf_fsav;
0199     end
0200     %
0201 end
0202 %
0203 j_rel_k = 0.003732;
0204 p_rel_k = 0.0001256;
0205 eta_rel_k = 29.72;
0206 %
0207 j_nr_k = 0.07092;
0208 p_nr_k = 0.004294;
0209 eta_nr_k = 16.52;
0210 %
0211 eta_rel_50 = j_rel_50./p_rel_50;
0212 eta_rel_100 = j_rel_100./p_rel_100;
0213 eta_rel_200 = j_rel_200./p_rel_200;
0214 eta_nr_50 = j_nr_50./p_nr_50;
0215 eta_nr_100 = j_nr_100./p_nr_100;
0216 eta_nr_200 = j_nr_200./p_nr_200;
0217 %
0218 %************************************************************************************************************************************
0219 %
0220 figure(1),clf
0221 %
0222 leg = {'n_{\xi} = 50','n_{\xi} = 100','n_{\xi} = 200'};
0223 xlim = 10.^[1,3];
0224 ylim = [0,0.01];
0225 xlab = 'n_p';
0226 ylab = 'j';
0227 tit = '';
0228 siz = 20+14i;
0229 %
0230 graph1D_jd(np_list,j_rel_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0231 graph1D_jd(np_list,j_rel_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0232 graph1D_jd(np_list,j_rel_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0233 graph1D_jd(xlim,[j_rel_k,j_rel_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0234 %
0235 set(gca,'ytick',[0:0.2:1]*ylim(2))
0236 set(gca,'xtick',[10,100,1000])
0237 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0238 set(gca,'XMinorGrid','off')
0239 set(gca,'XMinorTick','on')
0240 %
0241 figure(2),clf
0242 %
0243 ylim = [0,0.0005];
0244 ylab = 'P';
0245 %
0246 graph1D_jd(np_list,p_rel_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0247 graph1D_jd(np_list,p_rel_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0248 graph1D_jd(np_list,p_rel_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0249 graph1D_jd(xlim,[p_rel_k,p_rel_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0250 %
0251 set(gca,'ytick',[0:0.2:1]*ylim(2))
0252 set(gca,'xtick',[10,100,1000])
0253 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0254 set(gca,'XMinorGrid','off')
0255 set(gca,'XMinorTick','on')
0256 %
0257 figure(3),clf
0258 %
0259 ylim = [0,50];
0260 ylab = 'j/P';
0261 tit = '';
0262 %
0263 graph1D_jd(np_list,eta_rel_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0264 graph1D_jd(np_list,eta_rel_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0265 graph1D_jd(np_list,eta_rel_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0266 graph1D_jd(xlim,[eta_rel_k,eta_rel_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0267 %
0268 set(gca,'ytick',[0:0.2:1]*ylim(2))
0269 set(gca,'xtick',[10,100,1000])
0270 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0271 set(gca,'XMinorGrid','off')
0272 set(gca,'XMinorTick','on')
0273 %
0274 figure(4),clf
0275 %
0276 ylim = [0,0.1];
0277 ylab = 'j';
0278 tit = '';
0279 %
0280 graph1D_jd(np_list,j_nr_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0281 graph1D_jd(np_list,j_nr_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0282 graph1D_jd(np_list,j_nr_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0283 graph1D_jd(xlim,[j_nr_k,j_nr_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0284 %
0285 set(gca,'ytick',[0:0.2:1]*ylim(2))
0286 set(gca,'xtick',[10,100,1000])
0287 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0288 set(gca,'XMinorGrid','off')
0289 set(gca,'XMinorTick','on')
0290 %
0291 figure(5),clf
0292 %
0293 ylim = [0,0.01];
0294 ylab = 'P';
0295 %
0296 graph1D_jd(np_list,p_nr_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0297 graph1D_jd(np_list,p_nr_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0298 graph1D_jd(np_list,p_nr_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0299 graph1D_jd(xlim,[p_nr_k,p_nr_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0300 %
0301 set(gca,'ytick',[0:0.2:1]*ylim(2))
0302 set(gca,'xtick',[10,100,1000])
0303 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0304 set(gca,'XMinorGrid','off')
0305 set(gca,'XMinorTick','on')
0306 %
0307 figure(6),clf
0308 %
0309 ylim = [0,30];
0310 ylab = 'j/P';
0311 tit = '';
0312 %
0313 graph1D_jd(np_list,eta_nr_50,1,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca,0.9,0.7,0.7);
0314 graph1D_jd(np_list,eta_nr_100,1,0,'','','',NaN,xlim,ylim,'none','+','b',2,siz,gca);
0315 graph1D_jd(np_list,eta_nr_200,1,0,'','','',leg,xlim,ylim,'none','+','g',2,siz,gca);
0316 graph1D_jd(xlim,[eta_nr_k,eta_nr_k],0,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0317 %
0318 set(gca,'ytick',[0:0.2:1]*ylim(2))
0319 set(gca,'xtick',[10,100,1000])
0320 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0321 set(gca,'XMinorGrid','off')
0322 set(gca,'XMinorTick','on')
0323 %
0324 print_jd(p_opt,'fig_j_rel_grid','./figures',1)
0325 print_jd(p_opt,'fig_P_rel_grid','./figures',2)
0326 print_jd(p_opt,'fig_eta_rel_grid','./figures',3)
0327 print_jd(p_opt,'fig_j_nr_grid','./figures',4)
0328 print_jd(p_opt,'fig_P_nr_grid','./figures',5)
0329 print_jd(p_opt,'fig_eta_nr_grid','./figures',6)
0330 %
0331 %************************************************************************************************************************************
0332 %
0333 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0334 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.