Description of rundke

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_time';%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');dkeparam.invproc = -2;end
0061 %
0062 betath = 0.001;%validated for NR limit
0063 epsi = 0.04;%corresponds to pc/pT = 5
0064 %
0065 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0066 %
0067 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0068 equil.pzTi = ones(length(equil.zZi),1)*equil.pTe;
0069 %
0070 ohm = ohm_flat(equil,epsi);
0071 %
0072 RR_kulsrud = 1.914*1e-6;%Kulsrud (PRL, 31, 11, (1972) 690)
0073 %
0074 dkeparam.timevol = 1;%to calculate moments at all internal times
0075 %
0076 tn = logspace(0,5,101);
0077 nit = length(tn);
0078 dtn = diff([0,tn]);
0079 %
0080 dkeparam.dtn = NaN;% time steps prescribed by vector tn
0081 dkeparam.tn = tn;%when dtn is a vector of time steps
0082 %
0083 dkeparam.psin_S = psin_S;
0084 dkeparam.rho_S = rho_S;
0085 %
0086 dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0087 %
0088 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0089 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
0090 %
0091 [dummy,dummy,dummy,dke_out_0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0092 %
0093 xRRv_0 = NaN(1,nit);
0094 xRRs_0 = NaN(1,nit);
0095 for it=1:nit,
0096     xRRv_0(it) = (dke_out_0.normf0{it}(1) - dke_out_0.normf0{it}(end))/dtn(it);
0097     xRRs_0(it) = dke_out_0.XxRR_fsav{it}(end,:);
0098 end
0099 %
0100 dkeparam.boundary_mode_f = 1;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0101 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
0102 %
0103 [dummy,dummy,dummy,dke_out_1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0104 %
0105 xRRv_1 = NaN(1,nit);
0106 xRRs_1 = NaN(1,nit);
0107 for it=1:nit,
0108     xRRv_1(it) = (dke_out_1.normf0{it}(1) - dke_out_1.normf0{it}(end))/dtn(it);
0109     xRRs_1(it) = dke_out_1.XxRR_fsav{it}(end,:);
0110 end
0111 %
0112 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0113 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
0114 %
0115 [dummy,dummy,dummy,dke_out_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0116 %
0117 xRRv_2 = NaN(1,nit);
0118 xRRs_2 = NaN(1,nit);
0119 for it=1:nit,
0120     xRRv_2(it) = (dke_out_2.normf0{it}(1) - dke_out_2.normf0{it}(end))/dtn(it);
0121     xRRs_2(it) = dke_out_2.XxRR_fsav{it}(end,:);
0122 end
0123 %
0124 dkeparam.coll_mode = 1;% High-velocity limit
0125 %
0126 [dummy,dummy,dummy,dke_out_2_1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0127 %
0128 xRRv_2_1 = NaN(1,nit);
0129 xRRs_2_1 = NaN(1,nit);
0130 for it=1:nit,
0131     xRRv_2_1(it) = (dke_out_2_1.normf0{it}(1) - dke_out_2_1.normf0{it}(end))/dtn(it);
0132     xRRs_2_1(it) = dke_out_2_1.XxRR_fsav{it}(end,:);
0133 end
0134 %
0135 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0136 %
0137 [dummy,dummy,dummy,dke_out_2_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0138 %
0139 xresidue_2_2 = NaN(1,nit);
0140 xRRv_2_2 = NaN(1,nit);
0141 xRRs_2_2 = NaN(1,nit);
0142 for it=1:nit,
0143     xresidue_2_2(it) = dke_out_2_2.residu_f{it}(end);
0144     xRRv_2_2(it) = (dke_out_2_2.normf0{it}(1) - dke_out_2_2.normf0{it}(end))/dtn(it);
0145     xRRs_2_2(it) = dke_out_2_2.XxRR_fsav{it}(end,:);
0146 end
0147 xRRv_2_2(xresidue_2_2 > dkeparam.prec0_f) = NaN;
0148 xRRs_2_2(xresidue_2_2 > dkeparam.prec0_f) = NaN;
0149 %
0150 %************************************************************************************************************************************
0151 %
0152 figure(1),clf
0153 %
0154 leg = {'Free Mode','Fixed point','Normalization','Kulsrud'};
0155 xlim = 10.^[2,5];
0156 ylim = 10.^[-7,-5];
0157 xlab = 't/\tau_c';
0158 ylab = '\Gamma_R\tau_c';
0159 tit = '';
0160 siz = 20+14i;
0161 %
0162 graph1D_jd(tn,xRRs_0,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0163 graph1D_jd(tn,xRRs_1,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0164 graph1D_jd(tn,xRRs_2,1,1,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0165 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,1,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0166 %
0167 set(gca,'xtick',10.^[2:5])
0168 set(gca,'ytick',10.^[-7:-5])
0169 set(gca,'XMinorGrid','off')
0170 set(gca,'XMinorTick','on')
0171 set(gca,'YMinorGrid','off')
0172 set(gca,'YMinorTick','on')
0173 %
0174 figure(2),clf
0175 %
0176 leg = {'Linearized','High v limit','Maxwellian'};
0177 %
0178 graph1D_jd(tn,xRRs_2_2,1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0179 graph1D_jd(tn,xRRs_2_1,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0180 graph1D_jd(tn,xRRs_2,1,1,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0181 %
0182 set(gca,'xtick',10.^[2:5])
0183 set(gca,'ytick',10.^[-7:-5])
0184 set(gca,'XMinorGrid','off')
0185 set(gca,'XMinorTick','on')
0186 set(gca,'YMinorGrid','off')
0187 set(gca,'YMinorTick','on')
0188 %
0189 figure(3),clf
0190 xlim = 10.^[0,5];
0191 ylim = 10.^[-8,-4.99];
0192 %
0193 mask = xRRv_0 > 0;
0194 %
0195 leg = {'Volumic losses','Boundary losses','Runaway rate'};
0196 %
0197 graph1D_jd(tn(mask),xRRv_0(mask),1,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0198 graph1D_jd(tn,xRRs_0,1,1,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0199 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,1,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0200 graph1D_jd(tn(~mask),-xRRv_0(~mask),1,1,xlab,ylab,tit,NaN,xlim,ylim,'--','none','r',2,siz,gca);
0201 %
0202 set(gca,'ytick',[1e-010 1e-09 1e-08 1e-07 1e-06 1e-05])
0203 set(gca,'xtick',[1 10 100 1000 10000 100000])
0204 set(gca,'XMinorGrid','off')
0205 set(gca,'XMinorTick','on')
0206 set(gca,'YMinorGrid','off')
0207 set(gca,'YMinorTick','on')
0208 %
0209 print_jd(p_opt,'fig_runaway_time_vs','./figures',1)
0210 print_jd(p_opt,'fig_runaway_time_coll','./figures',2)
0211 print_jd(p_opt,'fig_runaway_time_mode','./figures',3)
0212 %
0213 %************************************************************************************************************************************
0214 %
0215 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0216 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
rundke_time

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
