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_betath';%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 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0063 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 0064 dkeparam.tn = 10000;%time for asymptotic solution with norm_mode_f = 1 0065 dkeparam.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn 0066 % 0067 dkeparam.psin_S = psin_S; 0068 dkeparam.rho_S = rho_S; 0069 % 0070 epsi = 0.04; 0071 % 0072 equil.pzTi = 1e-10*ones(size(equil.pzTi)); 0073 % 0074 betath_list = [logspace(-5,-2,7),logspace(-1.85,-0.35,11)]; 0075 % 0076 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants 0077 % 0078 nbetath = length(betath_list); 0079 % 0080 RR_0 = NaN(1,nbetath); 0081 RR_1 = NaN(1,nbetath); 0082 RR_2 = NaN(1,nbetath); 0083 %RR_4 = NaN(1,nbetath); 0084 % 0085 for ibetath = 1:nbetath, 0086 % 0087 betath = betath_list(ibetath); 0088 % 0089 equil.pTe = betath^2*mc2*ones(size(equil.pTe)); 0090 %equil.pzTi = betath^2*mc2*ones(size(equil.pzTi)); 0091 % 0092 ohm = ohm_flat(equil,epsi); 0093 % 0094 dkeparam.coll_mode = 0;% Relativistic Maxwellian background 0095 [dummy,dummy,dummy,dke_out_0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0096 % 0097 dkeparam.coll_mode = 1;% High-velocity limit 0098 [dummy,dummy,dummy,dke_out_1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0099 % 0100 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker 0101 [dummy,dummy,dummy,dke_out_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0102 % 0103 %dkeparam.coll_mode = 4;% Maxwellian Background 0104 %[dummy,dummy,dummy,dke_out_4] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0105 % 0106 RR_0(ibetath) = dke_out_0.XxRR_fsav(end,:); 0107 RR_1(ibetath) = dke_out_1.XxRR_fsav(end,:); 0108 % 0109 if dke_out_2.residu_f{end}(end) <= dkeparam.prec0_f, 0110 RR_2(ibetath) = dke_out_2.XxRR_fsav(end,:); 0111 end 0112 % 0113 % RR_4(ibetath) = dke_out_4.XxRR_fsav(end,:); 0114 % 0115 end 0116 % 0117 RR_kulsrud = 1.914*1e-6;%Kulsrud (PRL, 31,11, (1972) 690) 0118 % 0119 %************************************************************************************************************************************ 0120 % 0121 figure(1),clf 0122 % 0123 %leg = {'Linearized','High v limit','Maxwellian','Maxwellian NR'}; 0124 leg = {'Linearized','High v limit','Maxwellian'}; 0125 xlim = 10.^[-5,0]; 0126 ylim = [0,5e-6]; 0127 xlab = '\beta_T'; 0128 ylab = '\Gamma_R'; 0129 tit = ''; 0130 siz = 20+14i; 0131 % 0132 graph1D_jd(betath_list,RR_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7); 0133 graph1D_jd(betath_list,RR_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca); 0134 graph1D_jd(betath_list,RR_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca); 0135 %graph1D_jd(betath_list,RR_4,1,0,'','','',leg,xlim,ylim,'-','none','m',2,siz,gca); 0136 % 0137 graph1D_jd(xlim,[RR_kulsrud,RR_kulsrud],1,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca); 0138 set(gca,'ytick',[0:0.2:1]*ylim(2)) 0139 set(gca,'xtick',[1e-05 0.0001 0.001 0.01 0.1 1]) 0140 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'}) 0141 set(gca,'XMinorGrid','off') 0142 set(gca,'XMinorTick','on') 0143 % 0144 print_jd(p_opt,'fig_runaway_betath','./figures',1) 0145 % 0146 %************************************************************************************************************************************ 0147 % 0148 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']); 0149 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);