rundke_betath

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 %
0021 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0022 %
0023 id_simul = 'Ohm_betath';%Simulation ID
0024 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0025 %
0026 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0027 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)
0028 %
0029 id_path = '';%For all paths used by DKE solver
0030 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0031 %
0032 id_equil = 'TScyl';%For plasma equilibrium
0033 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0036 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_display = 'NO_DISPLAY';%For output code display
0039 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 id_ohm = '';%For Ohmic electric contribution
0042 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0045 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_transpfaste = '';%For fast electron radial transport
0048 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 id_ripple = '';%For fast electron magnetic ripple losses
0051 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0052 %
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %************************************************************************************************************************************
0056 %
0057 [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);
0058 %
0059 %************************************************************************************************************************************
0060 %
0061 if exist('dmumpsmex');dkeparam.invproc = -2;end
0062 %
0063 dkeparam.boundary_mode_f = 0;%Number of points where the Maxwellian distribution is enforced from p = 0 (p=0, free conservative mode but param_inv(1) must be less than 1e-4, otherwise 1e-3 is OK most of the time. Sensitive to the number of points in p)
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 = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0066 %
0067 dkeparam.psin_S = psin_S;
0068 dkeparam.rho_S = rho_S;
0069 %
0070 epsi = 0.001;
0071 %
0072 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0073 %
0074 betath_Karney = sqrt([0.01,0.02,0.05,0.1,0.2]);
0075 sigma_Karney = [7.291,7.160,6.807,6.317,5.575];
0076 sigma_Karney_nr_0 = 3.773;
0077 sigma_Karney_nr_1 = 2.837;
0078 sigma_Karney_nr_2 = 7.429;
0079 %sigma_Karney_nr_4 = 3.773;
0080 %
0081 betath_list = [logspace(-5,-2,7),0.015,0.02,0.03,0.05,0.07,betath_Karney];
0082 %
0083 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0084 %
0085 for ibetath = 1:length(betath_list),
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,Zcurr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0096     sigma_0(ibetath) = Zcurr.x_0/epsi;
0097     %
0098     dkeparam.coll_mode = 1;% High-velocity limit
0099     [dummy,Zcurr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0100     sigma_1(ibetath) = Zcurr.x_0/epsi;
0101     %
0102     dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0103     [dummy,Zcurr,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0104     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0105         sigma_2(ibetath) = Zcurr.x_0/epsi;
0106     else
0107         sigma_2(ibetath) = NaN;
0108     end        
0109     %
0110 end
0111 %
0112 sigma_Z1s = 1/100*sqrt(me*qe/10)/(1.65*mu0*e0^2);%Spitzer conductivity (from 'Tokamaks' by J. Wesson, Oxford Science Publication, 3rd Eds., 2004)
0113 %
0114 %************************************************************************************************************************************
0115 %
0116 format
0117 %
0118 delete res_betath
0119 %
0120 diary res_betath
0121 %
0122 disp(['Comparison LUKE/Karney/Spitzer - NR calculations (betath = 0.01)'])
0123 disp(['----------------------'])
0124 %
0125 disp(['Linearized operator : sigma = ',num2str(sigma_2(7)),'/',num2str(sigma_Karney_nr_2),'/',num2str(sigma_Z1s)]) 
0126 disp(['Maxwellian operator : sigma = ',num2str(sigma_0(7)),'/',num2str(sigma_Karney_nr_0)]) 
0127 disp(['High-Velocity operator : sigma = ',num2str(sigma_1(7)),'/',num2str(sigma_Karney_nr_1)]) 
0128 %
0129 disp(' ')
0130 disp(['Comparison LUKE/Karney - FR calculations'])
0131 disp(['----------------------'])
0132 disp(['betath^2 = 0.01 : sigma = ',num2str(sigma_2(13)),'/',num2str(sigma_Karney(1))]) 
0133 disp(['betath^2 = 0.02 : sigma = ',num2str(sigma_2(14)),'/',num2str(sigma_Karney(2))])
0134 disp(['betath^2 = 0.05 : sigma = ',num2str(sigma_2(15)),'/',num2str(sigma_Karney(3))])
0135 disp(['betath^2 = 0.10 : sigma = ',num2str(sigma_2(16)),'/',num2str(sigma_Karney(4))])
0136 disp(['betath^2 = 0.20 : sigma = ',num2str(sigma_2(17)),'/',num2str(sigma_Karney(5))])
0137 %
0138 diary off
0139 %
0140 figure(1),clf
0141 %
0142 %leg = {'Linearized','High v limit','Maxwellian NR','Maxwellian'};
0143 leg = {'Linearized','High v limit','Maxwellian'};
0144 xlim = 10.^[-5,0];
0145 ylim = [0,10];
0146 xlab = '\beta_T';
0147 ylab = '\sigma';
0148 tit = '';
0149 siz = 20+14i;
0150 %
0151 graph1D_jd(betath_list,sigma_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0152 graph1D_jd(betath_list,sigma_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0153 %graph1D_jd(betath_list,sigma_4,1,0,'','','',NaN,xlim,ylim,'-','none','m',2,siz,gca);
0154 graph1D_jd(betath_list,sigma_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0155 graph1D_jd(xlim,[sigma_Z1s,sigma_Z1s],1,0,'','','',NaN,xlim,ylim,'-','none','k',0.5,siz,gca);
0156 graph1D_jd(betath_Karney,sigma_Karney,0,0,xlab,ylab,tit,NaN,xlim,ylim,'none','+','r',2,siz,gca);
0157 graph1D_jd(xlim,[sigma_Karney_nr_2,sigma_Karney_nr_2],0,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0158 graph1D_jd(xlim,[sigma_Karney_nr_1,sigma_Karney_nr_1],0,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0159 %graph1D_jd(xlim,[sigma_Karney_nr_4,sigma_Karney_nr_4],0,0,'','','',NaN,xlim,ylim,'--','none','m',2,siz,gca);
0160 graph1D_jd(xlim,[sigma_Karney_nr_0,sigma_Karney_nr_0],0,0,'','','',NaN,xlim,ylim,'--','none','g',2,siz,gca);
0161 %
0162 set(gca,'ytick',[0:2:10])
0163 set(gca,'xtick',[1e-05 0.0001 0.001 0.01 0.1 1])
0164 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0165 set(gca,'XMinorGrid','off')
0166 set(gca,'XMinorTick','on')
0167 %
0168 print_jd(p_opt,'fig_sigma_betath','./figures',1)
0169 %
0170 %************************************************************************************************************************************
0171 %
0172 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0173 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);