rundke_time

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 = 'Ohm_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 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)
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.timevol = 1;%to calculate moments at all internal times
0065 %
0066 dkeparam.psin_S = psin_S;
0067 dkeparam.rho_S = rho_S;
0068 %
0069 nit = 141;
0070 tn = logspace(-2,5,nit);
0071 %
0072 dkeparam.tn = tn;
0073 %
0074 betath = 0.01;%validated for NR limit
0075 %
0076 epsi = 0.001;
0077 %
0078 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0079 %
0080 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0081 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0082 %
0083 ohm = ohm_flat(equil,epsi);
0084 %
0085 sigma_Karney_nr_0 = 3.773;
0086 sigma_Karney_nr_1 = 2.837;
0087 sigma_Karney_nr_2 = 7.429;
0088 %
0089 dkeparam.coll_mode = 0;% Relativistic Maxwellian background
0090 [dummy,Zcurr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0091 %
0092 sigma_0 = NaN(1,nit);
0093 for it=1:nit,
0094     sigma_0(it) = Zcurr(it).x_0/epsi;
0095 end
0096 %
0097 dkeparam.coll_mode = 1;% High-velocity limit
0098 [dummy,Zcurr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0099 %
0100 sigma_1 = NaN(1,nit);
0101 for it=1:nit,
0102     sigma_1(it) = Zcurr(it).x_0/epsi;
0103 end
0104 %
0105 dkeparam.coll_mode = 2;% Linearized Belaiev-Budker
0106 [dummy,Zcurr,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0107 %
0108 residue_2 = NaN(1,nit);
0109 sigma_2 = NaN(1,nit);
0110 for it=1:nit,
0111     residue_2(it) = dke_out.residu_f{it}(end);
0112     sigma_2(it) = Zcurr(it).x_0/epsi;
0113 end
0114 sigma_2(residue_2 > dkeparam.prec0_f) = NaN;
0115 %
0116 dkeparam.boundary_mode_f = 0;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0117 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
0118 [dummy,Zcurr,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0119 %
0120 residue_2_0 = NaN(1,nit);
0121 sigma_2_0 = NaN(1,nit);
0122 for it=1:nit,
0123     residue_2_0(it) = dke_out.residu_f{it}(end);
0124     sigma_2_0(it) = Zcurr(it).x_0/epsi;
0125 end
0126 sigma_2_0(residue_2_0 > dkeparam.prec0_f) = NaN;
0127 %
0128 dkeparam.boundary_mode_f = 1;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0129 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
0130 [dummy,Zcurr,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0131 %
0132 residue_2_1 = NaN(1,nit);
0133 sigma_2_1 = NaN(1,nit);
0134 for it=1:nit,
0135     residue_2_1(it) = dke_out.residu_f{it}(end);
0136     sigma_2_1(it) = Zcurr(it).x_0/epsi;
0137 end
0138 sigma_2_1(residue_2_1 > dkeparam.prec0_f) = NaN;
0139 %
0140 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)
0141 %
0142 %************************************************************************************************************************************
0143 %
0144 figure(1),clf
0145 %
0146 leg = {'Linearized','High v limit','Maxwellian'};
0147 xlim = 10.^[-2,5];
0148 ylim = [0,10];
0149 xlab = 't/\tau_c';
0150 ylab = '\sigma';
0151 tit = '';
0152 siz = 20+14i;
0153 %
0154 graph1D_jd(tn,sigma_2,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0155 graph1D_jd(tn,sigma_1,1,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0156 graph1D_jd(tn,sigma_0,1,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0157 graph1D_jd(xlim,[sigma_Karney_nr_2,sigma_Karney_nr_2],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0158 graph1D_jd(xlim,[sigma_Karney_nr_1,sigma_Karney_nr_1],1,0,'','','',NaN,xlim,ylim,'--','none','b',2,siz,gca);
0159 graph1D_jd(xlim,[sigma_Karney_nr_0,sigma_Karney_nr_0],1,0,'','','',NaN,xlim,ylim,'--','none','g',2,siz,gca);
0160 graph1D_jd(xlim,[sigma_Z1s,sigma_Z1s],1,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0161 %
0162 set(gca,'ytick',[0:2:10])
0163 set(gca,'xtick',10.^[-2:5])
0164 set(gca,'XMinorGrid','off')
0165 set(gca,'XMinorTick','on')
0166 %
0167 figure(2),clf
0168 %
0169 leg = {'Free Mode','Fixed point','Normalization'};
0170 %
0171 graph1D_jd(tn,sigma_2_0,1,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','b',2,siz,gca,0.9,0.7,0.7);
0172 graph1D_jd(tn,sigma_2_1,1,0,'','','',NaN,xlim,ylim,'-','none','g',2,siz,gca);
0173 graph1D_jd(tn,sigma_2,1,0,'','','',leg,xlim,ylim,'-','none','r',2,siz,gca);
0174 graph1D_jd(xlim,[sigma_Karney_nr_2,sigma_Karney_nr_2],1,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0175 graph1D_jd(xlim,[sigma_Z1s,sigma_Z1s],1,0,'','','',NaN,xlim,ylim,'--','none','k',2,siz,gca);
0176 %
0177 set(gca,'ytick',[0:2:10])
0178 set(gca,'xtick',10.^[-2:5])
0179 set(gca,'XMinorGrid','off')
0180 set(gca,'XMinorTick','on')
0181 %
0182 print_jd(p_opt,'fig_sigma_time','./figures',1)
0183 print_jd(p_opt,'fig_sigma_norm','./figures',2)
0184 %
0185 %************************************************************************************************************************************
0186 %
0187 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0188 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);