rundke_init

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 = 'Runaway_init';%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 if exist('dmumpsmex','file');dkeparam.invproc = -2;end
0061 %
0062 betath = 0.001;%validated for NR limit
0063 %
0064 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0065 %
0066 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0067 equil.pzTi = 1e-10*ones(size(equil.pzTi));
0068 %
0069 dkeparam.boundary_mode_f = 0;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0070 dkeparam.timevol = 1;%to calculate moments at all internal times
0071 %
0072 dkeparam.psin_S = psin_S;
0073 dkeparam.rho_S = rho_S;
0074 %
0075 dtn = 100;% time step
0076 nit = 100;% number of iterations
0077 tn_max = dtn*nit;% total time
0078 %
0079 epsi = 0.04;%corresponds to pc/pT = 5
0080 %
0081 ohm = ohm_flat(equil,epsi);
0082 %
0083 RR_kulsrud = 1.914*1e-6;%Kulsrud (PRL, 31,11, (1972) 690)
0084 %
0085 % ------------------------------
0086 % Internal time loop calculation
0087 % ------------------------------
0088 %
0089 dkeparam.dtn = dtn;    
0090 dkeparam.tn = tn_max;    
0091 %
0092 dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0093 %
0094 [Znorm_0,Zcurr_0,ZP0_0,dke_out_0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0095 %
0096 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0097 %
0098 [Znorm_2,Zcurr_2,ZP0_2,dke_out_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0099 %
0100 % One can check that dke_out_0.normf0{it}(1) is the same as dke_out_0.normf0{it-1}(end)
0101 %
0102 snormf0_0 = NaN(1,nit);
0103 sRRv_0 = NaN(1,nit);
0104 sRRs_0 = NaN(1,nit);
0105 for it=1:nit,
0106     snormf0_0(it) = dke_out_0.normf0{it}(end);
0107     sRRv_0(it) = (dke_out_0.normf0{it}(1) - dke_out_0.normf0{it}(end))/dtn;
0108     sRRs_0(it) = dke_out_0.XxRR_fsav{it}(end,:);
0109 end
0110 %
0111 snormf0_2 = NaN(1,nit);
0112 sresidue_2 = NaN(1,nit);
0113 sRRv_2 = NaN(1,nit);
0114 sRRs_2 = NaN(1,nit);
0115 for it=1:nit,
0116     snormf0_2(it) = dke_out_2.normf0{it}(end);
0117     sresidue_2(it) = dke_out_2.residu_f{it}(end);
0118     sRRv_2(it) = (dke_out_2.normf0{it}(1) - dke_out_2.normf0{it}(end))/dtn;
0119     sRRs_2(it) = dke_out_2.XxRR_fsav{it}(end,:);
0120 end
0121 snormf0_2(sresidue_2 > dkeparam.prec0_f) = NaN;
0122 sRRv_2(sresidue_2 > dkeparam.prec0_f) = NaN;
0123 sRRs_2(sresidue_2 > dkeparam.prec0_f) = NaN;
0124 %
0125 % ------------------------------
0126 % External time loop calculation
0127 % ------------------------------
0128 %
0129 XXf0_0 = [];
0130 XXf0_2 = [];
0131 %
0132 dkeparam.tn = dtn;%Enforcing a single iteration
0133 %
0134 xnormf0_0 = NaN(1,nit);
0135 xRRv_0 = NaN(1,nit);
0136 xRRs_0 = NaN(1,nit);
0137 %
0138 xnormf0_2 = NaN(1,nit);
0139 xRRv_2 = NaN(1,nit);
0140 xRRs_2 = NaN(1,nit);
0141 %
0142 for it = 1:nit,
0143     %
0144     dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0145     %
0146     [Znorm_0,Zcurr_0,ZP0_0,dke_out_0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_0,[]);
0147     %
0148     dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0149     %
0150     [Znorm_2,Zcurr_2,ZP0_2,dke_out_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,XXf0_2,[]);
0151     %
0152     % One can check that :
0153     %
0154     % - Znorm_0.x_0_fsav is the same as dke_out_0.normf0{1}(end)
0155     % - dke_out_0.normf0{1}(1) is the same as xnormf0_0(it-1)
0156     %
0157     xnormf0_0(it) = Znorm_0.x_0_fsav;
0158     xRRv_0(it) = (dke_out_0.normf0{1}(1) - dke_out_0.normf0{1}(end))/dtn;
0159     xRRs_0(it) = dke_out_0.XxRR_fsav{1}(end,:);
0160     %
0161     if dke_out_2.residu_f{1}(end) <= dkeparam.prec0_f,
0162         xnormf0_2(it) = Znorm_2.x_0_fsav;
0163         xRRv_2(it) = (dke_out_2.normf0{1}(1) - dke_out_2.normf0{1}(end))/dtn;
0164         xRRs_2(it) = dke_out_2.XxRR_fsav{1}(end,:);
0165     else
0166         xnormf0_2(it) = NaN;
0167         RRv_end_2 = NaN;
0168         RRs_end_2 = NaN;
0169     end        
0170     %
0171     XXf0_0 = dke_out_0.XXf0{1};
0172     XXf0_2 = dke_out_2.XXf0{1};
0173     %
0174 end
0175 %
0176 %************************************************************************************************************************************
0177 %
0178 format
0179 %
0180 delete res_init
0181 %
0182 diary res_init
0183 %
0184 disp('Test for internal and external loop consistency')
0185 disp('----------------------')
0186 disp('--> Maximum difference between internal and external density calculations : ')
0187 disp(['coll_mode = 0 : max(abs(snormf0_0 - xnormf0_0)) = ',num2str(max(abs(snormf0_0 - xnormf0_0)))]) 
0188 disp(['coll_mode = 2 : max(abs(snormf0_2 - xnormf0_2)) = ',num2str(max(abs(snormf0_2 - xnormf0_2)))]) 
0189 disp(' ')
0190 disp('--> Maximum difference between internal and external volume runaway calculations : ')
0191 disp(['coll_mode = 0 : max(abs(sRRv_0 - xRRv_0)) = ',num2str(max(abs(sRRv_0 - xRRv_0)))]) 
0192 disp(['coll_mode = 2 : max(abs(sRRv_2 - xRRv_2)) = ',num2str(max(abs(sRRv_2 - xRRv_2)))]) 
0193 disp(' ')
0194 disp('--> Maximum difference between internal and external surface runaway calculations : ')
0195 disp(['coll_mode = 0 : max(abs(sRRs_0 - xRRs_0)) = ',num2str(max(abs(sRRs_0 - xRRs_0)))]) 
0196 disp(['coll_mode = 2 : max(abs(sRRs_2 - xRRs_2)) = ',num2str(max(abs(sRRs_2 - xRRs_2)))]) 
0197 disp('----------------------')
0198 %
0199 diary off
0200 %
0201 figure(1),clf
0202 %
0203 leg = {'Single run','Multiple runs'};
0204 xlim = [0,tn_max];
0205 ylim = 10.^[-10,0];
0206 xlab = '\Deltat';
0207 ylab = '\Deltan/n_0';
0208 tit = '';
0209 siz = 20+14i;
0210 %
0211 graph1D_jd((1:nit)*dtn,1-snormf0_0,0,1,xlab,ylab,tit,NaN,xlim,ylim,'--','none','r',2,siz,gca,0.9,0.7,0.7);
0212 graph1D_jd((1:nit)*dtn,1-xnormf0_0,0,1,'','','',leg,xlim,ylim,'-','none','b',0.5,siz,gca);
0213 %
0214 set(gca,'ytick',[1e-010 1e-08 1e-06 1e-04 1e-02 1e-00])
0215 set(gca,'xtick',[0:tn_max/5:tn_max])
0216 set(gca,'XMinorGrid','off')
0217 %set(gca,'XMinorTick','on')
0218 %
0219 figure(2),clf
0220 %
0221 graph1D_jd((1:nit)*dtn,1-snormf0_2,0,1,xlab,ylab,tit,NaN,xlim,ylim,'--','none','r',2,siz,gca,0.9,0.7,0.7);
0222 graph1D_jd((1:nit)*dtn,1-xnormf0_2,0,1,'','','',leg,xlim,ylim,'-','none','b',0.5,siz,gca);
0223 %
0224 set(gca,'ytick',[1e-010 1e-08 1e-06 1e-04 1e-02 1e-00])
0225 set(gca,'xtick',[0:tn_max/5:tn_max])
0226 set(gca,'XMinorGrid','off')
0227 %set(gca,'XMinorTick','on')
0228 %
0229 figure(3),clf
0230 %
0231 leg = {'Volumic losses','Boundary losses','Runaway losses'};
0232 xlim = [dtn,tn_max];
0233 ylim = 10.^[-10,-5];
0234 xlab = '\Deltat';
0235 ylab = '\Gamma_R';
0236 tit = '';
0237 %
0238 graph1D_jd((1:nit)*dtn,sRRv_0,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0239 graph1D_jd((1:nit)*dtn,sRRs_0,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0240 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,1,'','','',leg,xlim,ylim,'--','none','k',2,20,gca);
0241 %
0242 set(gca,'ytick',[1e-010 1e-09 1e-08 1e-07 1e-06 1e-05])
0243 set(gca,'xtick',[1 10 100 1000 10000])
0244 set(gca,'XMinorGrid','off')
0245 set(gca,'XMinorTick','on')
0246 set(gca,'YMinorGrid','off')
0247 set(gca,'YMinorTick','on')
0248 %
0249 figure(4),clf
0250 %
0251 leg = {'Volumic losses','Boundary losses'};
0252 %
0253 graph1D_jd((1:nit)*dtn,sRRv_2,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0254 graph1D_jd((1:nit)*dtn,sRRs_2,1,1,'','','',leg,xlim,ylim,'-','none','b',2,siz,gca);
0255 %
0256 set(gca,'ytick',[1e-010 1e-09 1e-08 1e-07 1e-06 1e-05])
0257 set(gca,'xtick',[1 10 100 1000 10000])
0258 set(gca,'XMinorGrid','off')
0259 set(gca,'XMinorTick','on')
0260 set(gca,'YMinorGrid','off')
0261 set(gca,'YMinorTick','on')
0262 %
0263 print_jd(p_opt,'fig_runaway_init_cons_coll0','./figures',1)
0264 print_jd(p_opt,'fig_runaway_init_cons_coll2','./figures',2)
0265 print_jd(p_opt,'fig_runaway_init_RR_coll0','./figures',3)
0266 print_jd(p_opt,'fig_runaway_init_RR_coll2','./figures',4)
0267 %
0268 %************************************************************************************************************************************
0269 %
0270 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0271 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);