rundke_epsi_Ec

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 permission = test_permissions_yp;
0012 %
0013 if ~permission 
0014     disp('Please move the script to a local folder where you have write permission before to run it')
0015     return;
0016 end
0017 %
0018 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0019 %
0020 id_simul = 'Ohm_epsi';%Simulation ID
0021 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0022 %
0023 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0024 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)
0025 %
0026 id_path = '';%For all paths used by DKE solver
0027 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0028 %
0029 id_equil = 'TScyl';%For plasma equilibrium
0030 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0031 %
0032 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0033 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_display = 'NO_DISPLAY';%For output code display
0036 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_ohm = '';%For Ohmic electric contribution
0039 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0042 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 id_transpfaste = '';%For fast electron radial transport
0045 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_ripple = '';%For fast electron magnetic ripple losses
0048 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 %************************************************************************************************************************************
0051 %************************************************************************************************************************************
0052 %************************************************************************************************************************************
0053 %
0054 [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);
0055 %
0056 %************************************************************************************************************************************
0057 %
0058 if exist('dmumpsmex');dkeparam.invproc = -2;end
0059 %
0060 dkeparam.boundary_mode_f = 0;%2;%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)
0061 dkeparam.norm_mode_f = 0;%1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0062 %
0063 dkeparam.psin_S = psin_S;
0064 dkeparam.rho_S = rho_S;
0065 %
0066 betath = 0.1;%validated for NR limit
0067 %
0068 epsi_list = logspace(-2,-1,21);%logspace(-5,-1.4,19);
0069 %
0070 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0071 %
0072 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0073 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0074 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0075 %
0076 for iepsi = 1:length(epsi_list),
0077     %
0078     epsi = epsi_list(iepsi);
0079     %
0080     ohm = ohm_flat(equil,epsi);
0081     %
0082     dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0083     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0084     if dke_out.residu_f{end}(max([1,find(~isnan(dke_out.residu_f{end}))])) < dkeparam.prec0_f,
0085         sigma_0(iepsi) = Zcurr.x_0/epsi;
0086         sigma_0_noconv(iepsi) = NaN;
0087     else
0088         sigma_0(iepsi) = NaN;
0089         sigma_0_noconv(iepsi) = Zcurr.x_0/epsi;
0090     end        
0091     pn3 = momentumDKE.pn.^3;
0092     Xmhu = repmat(momentumDKE.mhu,[length(pn3),1]);
0093     dZcurrdpn = 2*pi*pn3.*integral_dke_jd(momentumDKE.dmhu,dke_out.XXf0(:,:,1).*Xmhu,2).'./Zmomcoef.gamma;
0094     ipnmax = find(dZcurrdpn == max(dZcurrdpn));
0095     if ipnmax < length(pn3),
0096         Ecmax_0(iepsi) = integral_dke_jd(momentumDKE.dpn,dZcurrdpn.*Zmomcoef.Ec_ref,2)./integral_dke_jd(momentumDKE.dpn,dZcurrdpn,2);
0097     else
0098         Ecmax_0(iepsi) = NaN;
0099     end
0100     %Ecmax_0(iepsi) = Zmomcoef.Ec_ref(ipnmax);
0101     %
0102     dkeparam.coll_mode = 1;% High-velocity limit
0103     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0104     if dke_out.residu_f{end}(max([1,find(~isnan(dke_out.residu_f{end}))])) < dkeparam.prec0_f,
0105         sigma_1(iepsi) = Zcurr.x_0/epsi;
0106         sigma_1_noconv(iepsi) = NaN;
0107     else
0108         sigma_1(iepsi) = NaN;
0109         sigma_1_noconv(iepsi) = Zcurr.x_0/epsi;
0110     end        
0111     pn3 = momentumDKE.pn.^3;
0112     Xmhu = repmat(momentumDKE.mhu,[length(pn3),1]);
0113     dZcurrdpn = 2*pi*pn3.*integral_dke_jd(momentumDKE.dmhu,dke_out.XXf0(:,:,1).*Xmhu,2).'./Zmomcoef.gamma;
0114     ipnmax = find(dZcurrdpn == max(dZcurrdpn));
0115     if ipnmax < length(pn3),
0116         Ecmax_1(iepsi) = integral_dke_jd(momentumDKE.dpn,dZcurrdpn.*Zmomcoef.Ec_ref,2)./integral_dke_jd(momentumDKE.dpn,dZcurrdpn,2);
0117     else
0118         Ecmax_1(iepsi) = NaN;
0119     end
0120     %Ecmax_1(iepsi) = Zmomcoef.Ec_ref(ipnmax);
0121     %
0122     dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0123     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0124     if dke_out.residu_f{end}(max([1,find(~isnan(dke_out.residu_f{end}))])) < dkeparam.prec0_f,
0125         sigma_2(iepsi) = Zcurr.x_0/epsi;
0126         sigma_2_noconv(iepsi) = NaN;
0127     else
0128         sigma_2(iepsi) = NaN;
0129         sigma_2_noconv(iepsi) = Zcurr.x_0/epsi;
0130     end        
0131     pn3 = momentumDKE.pn.^3;
0132     Xmhu = repmat(momentumDKE.mhu,[length(pn3),1]);
0133     dZcurrdpn = 2*pi*pn3.*integral_dke_jd(momentumDKE.dmhu,dke_out.XXf0(:,:,1).*Xmhu,2).'./Zmomcoef.gamma;
0134     ipnmax = find(dZcurrdpn == max(dZcurrdpn));
0135     if ipnmax < length(pn3),
0136         Ecmax_2(iepsi) = integral_dke_jd(momentumDKE.dpn,dZcurrdpn.*Zmomcoef.Ec_ref,2)./integral_dke_jd(momentumDKE.dpn,dZcurrdpn,2);
0137     else
0138         Ecmax_2(iepsi) = NaN;
0139     end
0140     %Ecmax_2(iepsi) = Zmomcoef.Ec_ref(ipnmax);
0141     %
0142 end
0143 %
0144 sigma_0_noconv(max(find(~isnan(sigma_0)))) = sigma_0(max(find(~isnan(sigma_0))));
0145 sigma_1_noconv(max(find(~isnan(sigma_1)))) = sigma_1(max(find(~isnan(sigma_1))));
0146 sigma_2_noconv(max(find(~isnan(sigma_2)))) = sigma_2(max(find(~isnan(sigma_2))));
0147 %
0148 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)
0149 %
0150 sigma_Karney_nr_0 = [3.773];
0151 sigma_Karney_nr_1 = [2.837];
0152 sigma_Karney_nr_2 = [7.429];
0153 %
0154 Ecmax_norm_0 = Ecmax_0;%/(mc2*betath.^2);
0155 Ecmax_norm_1 = Ecmax_1;%/(mc2*betath.^2);
0156 Ecmax_norm_2 = Ecmax_2;%/(mc2*betath.^2);
0157 %
0158 %************************************************************************************************************************************
0159 %
0160 figure(1),clf
0161 %
0162 leg = {'Linearized','High v limit','Maxwellian'};
0163 xlim = 10.^[-2,-1];%10.^[-5,0];%
0164 ylim = [0,10];
0165 xlab = 'E/E_D';
0166 ylab = '\sigma';
0167 tit = '';
0168 graph1D_jd(epsi_list,sigma_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,20,gca,0.9,0.7,0.7);
0169 graph1D_jd(epsi_list,sigma_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,20,gca);
0170 graph1D_jd(epsi_list,sigma_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,20,gca);
0171 graph1D_jd(epsi_list,sigma_2_noconv,1,0,'','','',NaN,xlim,ylim,'--','none','r',2,20,gca);
0172 graph1D_jd(epsi_list,sigma_1_noconv,1,0,'','','',NaN,xlim,ylim,'--','none','b',2,20,gca);
0173 graph1D_jd(epsi_list,sigma_0_noconv,1,0,'','','',NaN,xlim,ylim,'--','none','g',2,20,gca);
0174 graph1D_jd(xlim,[sigma_Z1s,sigma_Z1s],1,0,'','','',NaN,xlim,ylim,'-','none','k',0.5,20,gca);
0175 graph1D_jd(xlim,[sigma_Karney_nr_2,sigma_Karney_nr_2],0,0,'','','',NaN,xlim,ylim,'--','none','r',0.5,20,gca);
0176 graph1D_jd(xlim,[sigma_Karney_nr_1,sigma_Karney_nr_1],0,0,'','','',NaN,xlim,ylim,'--','none','b',0.5,20,gca);
0177 graph1D_jd(xlim,[sigma_Karney_nr_0,sigma_Karney_nr_0],0,0,'','','',NaN,xlim,ylim,'--','none','g',0.5,20,gca);
0178 graph1D_jd([1,1]/dkeparam.pnmax_S^2,ylim,0,0,'','','',NaN,xlim,ylim,'--','none','k',2,20,gca);
0179 set(gca,'ytick',[0:2:10])
0180 set(gca,'xtick',[1e-06 1e-5 0.0001 0.001 0.01 0.1 1])
0181 %set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0182 set(gca,'XMinorGrid','off')
0183 set(gca,'XMinorTick','on')
0184 %
0185 print_jd(2,'fig_sigma_epsi','./figures',1)
0186 %
0187 figure(2),clf
0188 %
0189 leg = {'Linearized','High v limit','Maxwellian'};
0190 ylim = [0,200];
0191 xlab = 'E/E_D';
0192 ylab = '<E> (keV)';%'E_{Jmax}/T';
0193 tit = '';
0194 graph1D_jd(epsi_list,Ecmax_norm_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,20,gca,0.9,0.7,0.7);
0195 %graph1D_jd(epsi_list,Ecmax_norm_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,20,gca);
0196 %graph1D_jd(epsi_list,Ecmax_norm_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,20,gca);
0197 set(gca,'ytick',[0:40:200])
0198 set(gca,'xtick',[0.01:0.01:0.1])%set(gca,'xtick',[1e-06 1e-5 0.0001 0.001 0.01 0.1 1])
0199 set(gca,'XTickLabel',['0.01',cell(1,8),'0.1'])%set(gca,'XTickLabel',{'10^{-5}','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})
0200 set(gca,'XMinorGrid','off')
0201 set(gca,'XMinorTick','on')
0202 %
0203 print_jd(2,'fig_Ecmax_epsi','./figures',2)
0204 %
0205 %************************************************************************************************************************************
0206 %
0207 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0208 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);