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 in $HOME/Database directory if necessary *****************************
0021 %
0022 id_simul = 'Thermal_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_display = 'NO_DISPLAY';%For output code display
0035 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0036 %
0037 id_ohm = '';%For Ohmic electric contribution
0038 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0039 %
0040 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0041 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0042 %
0043 id_transpfaste = '';%For fast electron radial transport
0044 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0045 %
0046 id_ripple = '';%For fast electron magnetic ripple losses
0047 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0048 %
0049 %************************************************************************************************************************************
0050 %************************************************************************************************************************************
0051 %************************************************************************************************************************************
0052 %
0053 [dkepath,equil,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'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);
0054 %
0055 %************************************************************************************************************************************
0056 %
0057 if exist('dmumpsmex','file'),
0058     invproc = -2;% MUMPSMEX
0059 else
0060     invproc = 1;
0061 end
0062 %
0063 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 coll_mode = 2;%Relativistic collision operator: (0) Relativistic Maxwellian background, (1) High-velocity limit, (2) Linearized Belaiev-Budker (3) Non-relativistic Lorentz model, (4) Non-relativistic Maxwellian background
0065 %
0066 tnmax = 100000;%asymptotic time
0067 dtn_list = [500,5000,50000];%time steps
0068 %
0069 pzTi = equil.pzTi;
0070 %
0071 % Case Ti=Te, uniform grid
0072 %
0073 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0074 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0075 %
0076 [dkeparam] = load_structures_yp('dkeparam',id_dkeparam,path_dkeparam);
0077 %
0078 dkeparam.boundary_mode_f = boundary_mode_f;
0079 dkeparam.coll_mode = coll_mode;
0080 dkeparam.invproc = invproc;
0081 %
0082 dkeparam.prec0_f = 0;%to reach end of Legendre iterations
0083 dkeparam.tn = tnmax;
0084 %
0085 dkeparam.psin_S = psin_S;
0086 dkeparam.rho_S = rho_S;
0087 %
0088 ndtn = length(dtn_list);
0089 %
0090 normf0_u_1 = cell(1,ndtn);
0091 normf0_u_2 = cell(1,ndtn);
0092 normf0_nu_1 = cell(1,ndtn);
0093 normf0_nu_2 = cell(1,ndtn);
0094 %
0095 residue_u_1 = cell(1,ndtn);
0096 residue_u_2 = cell(1,ndtn);
0097 residue_nu_1 = cell(1,ndtn);
0098 residue_nu_2 = cell(1,ndtn);
0099 %
0100 for idtn = 1:ndtn,
0101     %
0102     dtn = dtn_list(idtn);
0103     dkeparam.dtn = dtn;
0104     %
0105     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0106     %
0107     for itn = 1:length(dke_out.residu_f),
0108         normf0_u_1{idtn}(itn) = dke_out.normf0{itn}(end);
0109         residue_u_1{idtn}(itn) = dke_out.residu_f{itn}(end);
0110     end
0111 end
0112 %
0113 % Case Ti=Te/2, uniform grid
0114 %
0115 equil.pzTi = pzTi/2;
0116 %
0117 for idtn = 1:ndtn;
0118     %
0119     dtn = dtn_list(idtn);
0120     dkeparam.dtn = dtn;
0121     %
0122     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0123     %
0124     for itn = 1:length(dke_out.residu_f),
0125         normf0_u_2{idtn}(itn) = dke_out.normf0{itn}(end);
0126         residue_u_2{idtn}(itn) = dke_out.residu_f{itn}(end);
0127     end
0128 end
0129 %
0130 % Case Ti=Te, non-uniform grid
0131 %
0132 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0133 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0134 %
0135 [dkeparam] = load_structures_yp('dkeparam',id_dkeparam,path_dkeparam);
0136 %
0137 dkeparam.boundary_mode_f = boundary_mode_f;
0138 dkeparam.coll_mode = coll_mode;
0139 dkeparam.invproc = invproc;
0140 %
0141 dkeparam.prec0_f = 0;%to reach end of Legendre iterations
0142 dkeparam.tn = tnmax;
0143 %
0144 dkeparam.psin_S = psin_S;
0145 dkeparam.rho_S = rho_S;
0146 %
0147 equil.pzTi = pzTi;
0148 %
0149 for idtn = 1:ndtn;
0150     %
0151     dtn = dtn_list(idtn);
0152     dkeparam.dtn = dtn;
0153     %
0154     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0155     %
0156     for itn = 1:length(dke_out.residu_f),
0157         normf0_nu_1{idtn}(itn) = dke_out.normf0{itn}(end);
0158         residue_nu_1{idtn}(itn) = dke_out.residu_f{itn}(end);
0159     end
0160 end
0161 %
0162 % Case Ti=Te/2, non-uniform grid
0163 %
0164 equil.pzTi = pzTi/2;
0165 %
0166 for idtn = 1:ndtn;
0167     %
0168     dtn = dtn_list(idtn);
0169     dkeparam.dtn = dtn;
0170     %
0171     [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0172     %
0173     for itn = 1:length(dke_out.residu_f),
0174         normf0_nu_2{idtn}(itn) = dke_out.normf0{itn}(end);
0175         residue_nu_2{idtn}(itn) = dke_out.residu_f{itn}(end);
0176     end
0177 end
0178 %
0179 % density conservation
0180 %
0181 figure(1),clf
0182 %
0183 leg = {'\Deltat/\tau_c = 500','\Deltat/\tau_c = 5000','\Deltat/\tau_c = 50000','Location','SouthWest'};
0184 xlim = [0,tnmax];
0185 ylim = [-4,4]*1e-8;
0186 xlab = 't/\tau_c';
0187 ylab = '\Deltan/n_0';
0188 tit = '';
0189 colors = {'r','b','g'};
0190 siz = 20+14i;
0191 %
0192 for idtn = 1:ndtn-1;
0193     graph1D_jd(0:dtn_list(idtn):tnmax,[0,normf0_u_1{idtn}-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{idtn},2,siz,gca);
0194 end
0195 graph1D_jd(0:dtn_list(end):tnmax,[0,normf0_u_1{end}-1],0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{end},2,siz,gca,0.9,0.7,0.7);
0196 %
0197 set(gca,'xtick',[0:0.2:1]*tnmax)
0198 set(gca,'ytick',[-4:2:4]*1e-8)
0199 %
0200 figure(2),clf
0201 %
0202 ylim = [-4,4]*1e-3;
0203 %
0204 for idtn = 1:ndtn-1;
0205     graph1D_jd(0:dtn_list(idtn):tnmax,[0,normf0_u_2{idtn}-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{idtn},2,siz,gca);
0206 end
0207 graph1D_jd(0:dtn_list(end):tnmax,[0,normf0_u_2{end}-1],0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{end},2,siz,gca,0.9,0.7,0.7);
0208 %
0209 set(gca,'xtick',[0:0.2:1]*tnmax)
0210 set(gca,'ytick',[-4:2:4]*1e-3)
0211 %
0212 figure(3),clf
0213 %
0214 ylim = [-4,4]*1e-8;
0215 %
0216 for idtn = 1:ndtn-1;
0217     graph1D_jd(0:dtn_list(idtn):tnmax,[0,normf0_nu_1{idtn}-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{idtn},2,siz,gca);
0218 end
0219 graph1D_jd(0:dtn_list(end):tnmax,[0,normf0_nu_1{end}-1],0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{end},2,siz,gca,0.9,0.7,0.7);
0220 %
0221 set(gca,'xtick',[0:0.2:1]*tnmax)
0222 set(gca,'ytick',[-4:2:4]*1e-8)
0223 %
0224 figure(4),clf
0225 %
0226 ylim = [-4,4]*1e-3;
0227 %
0228 for idtn = 1:ndtn-1;
0229     graph1D_jd(0:dtn_list(idtn):tnmax,[0,normf0_nu_2{idtn}-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{idtn},2,siz,gca);
0230 end
0231 graph1D_jd(0:dtn_list(end):tnmax,[0,normf0_nu_2{end}-1],0,0,xlab,ylab,tit,leg,xlim,ylim,'-','none',colors{end},2,siz,gca,0.9,0.7,0.7);
0232 %
0233 set(gca,'xtick',[0:0.2:1]*tnmax)
0234 set(gca,'ytick',[-4:2:4]*1e-3)
0235 %
0236 % residue
0237 %
0238 figure(5),clf
0239 %
0240 ylim = [1e-24,2e-21];
0241 ytick = 10.^[-24:-21];
0242 ylab = 'Legendre convergence parameter';
0243 %
0244 for idtn = 1:ndtn-1;
0245     graph1D_jd(dtn_list(idtn):dtn_list(idtn):tnmax,residue_u_1{idtn},0,1,'','','',NaN,xlim,ylim,'-','+',colors{idtn},2,siz,gca);
0246 end
0247 graph1D_jd(dtn_list(end):dtn_list(end):tnmax,residue_u_1{end},0,1,xlab,ylab,tit,leg,xlim,ylim,'-','+',colors{end},2,siz,gca,0.9,0.7,0.7);
0248 %
0249 set(gca,'xtick',[0:0.2:1]*tnmax)
0250 set(gca,'ytick',ytick)
0251 %
0252 figure(6),clf
0253 %
0254 for idtn = 1:ndtn-1;
0255     graph1D_jd(dtn_list(idtn):dtn_list(idtn):tnmax,residue_u_2{idtn},0,1,'','','',NaN,xlim,ylim,'-','+',colors{idtn},2,siz,gca);
0256 end
0257 graph1D_jd(dtn_list(end):dtn_list(end):tnmax,residue_u_2{end},0,1,xlab,ylab,tit,leg,xlim,ylim,'-','+',colors{end},2,siz,gca,0.9,0.7,0.7);
0258 %
0259 set(gca,'xtick',[0:0.2:1]*tnmax)
0260 set(gca,'ytick',ytick)
0261 %
0262 figure(7),clf
0263 %
0264 for idtn = 1:ndtn-1;
0265     graph1D_jd(dtn_list(idtn):dtn_list(idtn):tnmax,residue_nu_1{idtn},0,1,'','','',NaN,xlim,ylim,'-','+',colors{idtn},2,siz,gca);
0266 end
0267 graph1D_jd(dtn_list(end):dtn_list(end):tnmax,residue_nu_1{end},0,1,xlab,ylab,tit,leg,xlim,ylim,'-','+',colors{end},2,siz,gca,0.9,0.7,0.7);
0268 %
0269 set(gca,'xtick',[0:0.2:1]*tnmax)
0270 set(gca,'ytick',ytick)
0271 %
0272 figure(8),clf
0273 %
0274 for idtn = 1:ndtn-1;
0275     graph1D_jd(dtn_list(idtn):dtn_list(idtn):tnmax,residue_nu_2{idtn},0,1,'','','',NaN,xlim,ylim,'-','+',colors{idtn},2,siz,gca);
0276 end
0277 graph1D_jd(dtn_list(end):dtn_list(end):tnmax,residue_nu_2{end},0,1,xlab,ylab,tit,leg,xlim,ylim,'-','+',colors{end},2,siz,gca,0.9,0.7,0.7);
0278 %
0279 set(gca,'xtick',[0:0.2:1]*tnmax)
0280 set(gca,'ytick',ytick)
0281 %
0282 print_jd(p_opt,'fig_thermal_dtn_TiTe_u','./figures',1)
0283 print_jd(p_opt,'fig_thermal_dtn_TiTe2_u','./figures',2)
0284 print_jd(p_opt,'fig_thermal_dtn_TiTe_nu','./figures',3)
0285 print_jd(p_opt,'fig_thermal_dtn_TiTe2_nu','./figures',4)
0286 %
0287 print_jd(p_opt,'fig_thermal_dtn_TiTe_u_res','./figures',5)
0288 print_jd(p_opt,'fig_thermal_dtn_TiTe2_u_res','./figures',6)
0289 print_jd(p_opt,'fig_thermal_dtn_TiTe_nu_res','./figures',7)
0290 print_jd(p_opt,'fig_thermal_dtn_TiTe2_nu_res','./figures',8)
0291 %
0292 %************************************************************************************************************************************
0293 %
0294 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0295 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);