rundke_Zi

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 = 'Runaway_Zi_nobounce';%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 %
0061 if exist('dmumpsmex');dkeparam.invproc = -2;end
0062 %
0063 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0064 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
0065 dkeparam.tn = 10000;%time for asymptotic solution with norm_mode_f = 1
0066 dkeparam.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn
0067 %
0068 dkeparam.psin_S = psin_S;
0069 dkeparam.rho_S = rho_S;
0070 %
0071 dkeparam.coll_mode = 0;%For comparison with theoretical results
0072 betath = 0.001;%validated for NR limit
0073 %
0074 equil.pzTi = 1e-10*ones(size(equil.pzTi));
0075 %
0076 Zi_list = [1:20];
0077 epsi_list = [0.04,0.06,0.08,0.10];
0078 %
0079 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0080 %
0081 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0082 %
0083 for iZi = 1:length(Zi_list),
0084     %
0085     Zi = Zi_list(iZi);
0086     equil.zZi = [1,1,1,Zi];
0087     equil.pzni = [0;0;0;1/Zi]*equil.pne;
0088     %
0089     for iepsi = 1:length(epsi_list),
0090         %
0091         epsi = epsi_list(iepsi);
0092         %
0093         ohm = ohm_flat(equil,epsi);
0094         %
0095         [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0096         %
0097         RR(iZi,iepsi) = dke_out.XxRR_fsav(end,:);
0098         %
0099     end
0100     %
0101 end
0102 %
0103 Zi_kulsrud =[1,2,3,10];
0104 RR_kulsrud = 1e-4*[ 0.01914,NaN   ,NaN  ,NaN  ;...%Kulsrud (PRL, 31,11, (1972) 690)
0105                     0.5411 ,0.2611,NaN  ,NaN  ;...  
0106                     3.177  ,1.735 ,1.047,0.09 ;...
0107                     10.04  ,5.839 ,3.757,0.449].';
0108 %
0109 %************************************************************************************************************************************
0110 %
0111 format
0112 %
0113 delete res_Zi
0114 %
0115 diary res_Zi
0116 %
0117 disp(['Comparison LUKE/Kulsrud - NR calculations (betath = 0.001)'])
0118 disp(['----------------------'])
0119 disp(['--> epsi = 0.04'])
0120 disp(['Zi =  1 : gamma_R = ',num2str(RR( 1,1)),'/',num2str(RR_kulsrud(1,1))]) 
0121 disp(['Zi =  2 : gamma_R = ',num2str(RR( 2,1)),'/',num2str(RR_kulsrud(2,1))]) 
0122 disp(['Zi =  3 : gamma_R = ',num2str(RR( 3,1)),'/',num2str(RR_kulsrud(3,1))]) 
0123 disp(['Zi = 10 : gamma_R = ',num2str(RR(10,1)),'/',num2str(RR_kulsrud(4,1))]) 
0124 disp(['----------------------'])
0125 disp(['--> epsi = 0.06'])
0126 disp(['Zi =  1 : gamma_R = ',num2str(RR( 1,2)),'/',num2str(RR_kulsrud(1,2))]) 
0127 disp(['Zi =  2 : gamma_R = ',num2str(RR( 2,2)),'/',num2str(RR_kulsrud(2,2))]) 
0128 disp(['Zi =  3 : gamma_R = ',num2str(RR( 3,2)),'/',num2str(RR_kulsrud(3,2))]) 
0129 disp(['Zi = 10 : gamma_R = ',num2str(RR(10,2)),'/',num2str(RR_kulsrud(4,2))]) 
0130 disp(['----------------------'])
0131 disp(['--> epsi = 0.08'])
0132 disp(['Zi =  1 : gamma_R = ',num2str(RR( 1,3)),'/',num2str(RR_kulsrud(1,3))]) 
0133 disp(['Zi =  2 : gamma_R = ',num2str(RR( 2,3)),'/',num2str(RR_kulsrud(2,3))]) 
0134 disp(['Zi =  3 : gamma_R = ',num2str(RR( 3,3)),'/',num2str(RR_kulsrud(3,3))]) 
0135 disp(['Zi = 10 : gamma_R = ',num2str(RR(10,3)),'/',num2str(RR_kulsrud(4,3))]) 
0136 disp(['----------------------'])
0137 disp(['--> epsi = 0.10'])
0138 disp(['Zi =  1 : gamma_R = ',num2str(RR( 1,4)),'/',num2str(RR_kulsrud(1,4))]) 
0139 disp(['Zi =  2 : gamma_R = ',num2str(RR( 2,4)),'/',num2str(RR_kulsrud(2,4))]) 
0140 disp(['Zi =  3 : gamma_R = ',num2str(RR( 3,4)),'/',num2str(RR_kulsrud(3,4))]) 
0141 disp(['Zi = 10 : gamma_R = ',num2str(RR(10,4)),'/',num2str(RR_kulsrud(4,4))]) 
0142 %
0143 diary off
0144 %
0145 figure(1),clf
0146 %
0147 leg = {'E/E_D = 0.04','E/E_D = 0.06','E/E_D = 0.08','E/E_D = 0.10'};
0148 xlim = [0,20];
0149 ylim = 10.^[-8,-2];
0150 xlab = 'Z_i';
0151 ylab = '\Gamma_R';
0152 tit = '';
0153 siz = 20+14i;
0154 %
0155 graph1D_jd(Zi_list,RR,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','',2,siz,gca,0.9,0.7,0.7);
0156 graph1D_jd(Zi_kulsrud,RR_kulsrud,0,1,'','','',leg,xlim,ylim,'none','+','',2,siz,gca);
0157 %
0158 set(gca,'ytick',[1e-8 1e-7 1e-6 1e-5 1e-4 1e-3 1e-2])
0159 set(gca,'xtick',[0:2:20])
0160 set(gca,'YMinorGrid','off')
0161 set(gca,'YMinorTick','on')
0162 %
0163 print_jd(p_opt,'fig_runaway_Zi','./figures',1)
0164 %
0165 %************************************************************************************************************************************
0166 %
0167 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0168 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.