rundke_dtn

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_dtn';%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.coll_mode = 0;% Relativistic Maxwellian background
0070 %
0071 dkeparam.psin_S = psin_S;
0072 dkeparam.rho_S = rho_S;
0073 %
0074 epsi = 0.04;%corresponds to pc/pT = 5
0075 %
0076 ohm = ohm_flat(equil,epsi);
0077 %
0078 RR_kulsrud = 1.914*1e-6;%Kulsrud (PRL, 31,11, (1972) 690)
0079 %
0080 % Building the distribution function to a quasi steady-state for the
0081 % runaway problem : tn = 10000
0082 %
0083 tnmax = 10000;
0084 nit_list = round(logspace(0,2,11));
0085 dtn_list = tnmax./nit_list;
0086 %
0087 dkeparam.tn = NaN;%specified by dtn
0088 %
0089 % Testing different dtn
0090 %
0091 ndtn = length(dtn_list);
0092 xRRv_0 = NaN(1,ndtn);
0093 xRRs_0 = NaN(1,ndtn);
0094 xRRv_1 = NaN(1,ndtn);
0095 xRRs_1 = NaN(1,ndtn);
0096 xRRv_2 = NaN(1,ndtn);
0097 xRRs_2 = NaN(1,ndtn);
0098 %
0099 for idtn = 1:ndtn,
0100     %
0101     dtn = dtn_list(idtn); 
0102     dkeparam.dtn = repmat(dtn,[1,nit_list(idtn)]);
0103     %
0104     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
0105     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
0106     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0107     %
0108     xRRv_0(idtn) = (dke_out.normf0{end}(1) - dke_out.normf0{end}(end))/dtn;
0109     xRRs_0(idtn) = dke_out.XxRR_fsav(end,:);
0110     %
0111     %
0112     dkeparam.boundary_mode_f = 1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0113     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
0114     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0115     %
0116     xRRv_1(idtn) = (dke_out.normf0{end}(1) - dke_out.normf0{end}(end))/dtn;
0117     xRRs_1(idtn) = dke_out.XxRR_fsav(end,:);
0118     %
0119     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
0120     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
0121     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0122     %
0123     xRRv_2(idtn) = (dke_out.normf0{end}(1) - dke_out.normf0{end}(end))/dtn;
0124     xRRs_2(idtn) = dke_out.XxRR_fsav(end,:);
0125     %
0126 end
0127 %
0128 %************************************************************************************************************************************
0129 %
0130 figure(1),clf
0131 %
0132 leg = {'Volumic losses','Boundary losses','Runaway losses'};
0133 xlim = 10.^[2,4];
0134 ylim = [1.5,2]*1e-6;
0135 xlab = '\Deltat';
0136 ylab = '\Gamma_R';
0137 tit = '';
0138 siz = 20+14i;
0139 %
0140 graph1D_jd(dtn_list,xRRv_0,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0141 graph1D_jd(dtn_list,xRRs_0,1,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0142 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0143 %
0144 set(gca,'ytick',[1.5:0.1:2]*1e-6)
0145 %
0146 figure(2),clf
0147 %
0148 graph1D_jd(dtn_list,xRRv_1,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0149 graph1D_jd(dtn_list,xRRs_1,1,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0150 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0151 %
0152 set(gca,'ytick',[1.5:0.1:2]*1e-6)
0153 %
0154 figure(3),clf
0155 %
0156 graph1D_jd(dtn_list,xRRv_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0157 graph1D_jd(dtn_list,xRRs_2,1,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0158 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0159 %
0160 set(gca,'ytick',[1.5:0.1:2]*1e-6)
0161 %
0162 print_jd(p_opt,'fig_runaway_dtn_0','./figures',1)
0163 print_jd(p_opt,'fig_runaway_dtn_1','./figures',2)
0164 print_jd(p_opt,'fig_runaway_dtn_2','./figures',3)
0165 %
0166 %************************************************************************************************************************************
0167 %
0168 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0169 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);