rundke_norm

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 in $HOME/Database directory if necessary *****************************
0019 %
0020 id_simul = 'Thermal_norm';%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 dtn = 10000;
0061 nit = 10;
0062 tnmax = dtn*nit;
0063 %
0064 dkeparam.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 dkeparam.dtn = dtn;
0066 dkeparam.tn = tnmax;
0067 dkeparam.prec0_f = 0;%to reach end of Legendre iterations
0068 %
0069 % Case Ti=Te, free mode
0070 %
0071 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)
0072 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
0073 %
0074 dkeparam.psin_S = psin_S;
0075 dkeparam.rho_S = rho_S;
0076 %
0077 normf0_1_0 = NaN(1,nit);
0078 normf0_1_1 = NaN(1,nit);
0079 normf0_1_2 = NaN(1,nit);
0080 normf0_2_0 = NaN(1,nit);
0081 normf0_2_1 = NaN(1,nit);
0082 normf0_2_2 = NaN(1,nit);
0083 %
0084 residue_1_0 = NaN(1,nit);
0085 residue_1_1 = NaN(1,nit);
0086 residue_1_2 = NaN(1,nit);
0087 residue_2_0 = NaN(1,nit);
0088 residue_2_1 = NaN(1,nit);
0089 residue_2_2 = NaN(1,nit);
0090 %
0091 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0092 %
0093 for itn = 1:nit,
0094     normf0_1_0(itn) = dke_out.normf0{itn}(end);
0095     residue_1_0(itn) = dke_out.residu_f{itn}(end);
0096 end
0097 %
0098 % Case Ti=Te, boundary mode
0099 %
0100 dkeparam.boundary_mode_f = 1;%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)
0101 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
0102 %
0103 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0104 %
0105 for itn = 1:nit,
0106     normf0_1_1(itn) = dke_out.normf0{itn}(end);
0107     residue_1_1(itn) = dke_out.residu_f{itn}(end);
0108 end
0109 %
0110 % Case Ti=Te, norm mode
0111 %
0112 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)
0113 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
0114 %
0115 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0116 %
0117 for itn = 1:nit,
0118     normf0_1_2(itn) = dke_out.normf0{itn}(end);
0119     residue_1_2(itn) = dke_out.residu_f{itn}(end);
0120 end
0121 %
0122 % Case Ti=Te/2, free mode
0123 %
0124 equil.pzTi = equil.pzTi/2;
0125 %
0126 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)
0127 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
0128 %
0129 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0130 %
0131 for itn = 1:nit,
0132     normf0_2_0(itn) = dke_out.normf0{itn}(end);
0133     residue_2_0(itn) = dke_out.residu_f{itn}(end);
0134 end
0135 %
0136 % Case Ti=Te/2, boundary mode
0137 %
0138 dkeparam.boundary_mode_f = 1;%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)
0139 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
0140 %
0141 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0142 %
0143 for itn = 1:nit,
0144     normf0_2_1(itn) = dke_out.normf0{itn}(end);
0145     residue_2_1(itn) = dke_out.residu_f{itn}(end);
0146 end
0147 %
0148 % Case Ti=Te/2, norm mode
0149 %
0150 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)
0151 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
0152 %
0153 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0154 %
0155 for itn = 1:nit,
0156     normf0_2_2(itn) = dke_out.normf0{itn}(end);
0157     residue_2_2(itn) = dke_out.residu_f{itn}(end);
0158 end
0159 %
0160 figure(1),clf
0161 %
0162 leg = {'Free Mode (\times 10^3)','Fixed point','Normalization (\times 10^3)'};
0163 xlim = [0,tnmax];
0164 ylim = [-8,8]*1e-5;
0165 xlab = 't/\tau_c';
0166 ylab = '\Deltan/n_0';
0167 tit = '';
0168 colors = {'r','b','g','m'};
0169 siz = 20+14i;
0170 %
0171 graph1D_jd(0:dtn:tnmax,[0,(normf0_1_0-1)*1e3],0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none',colors{1},2,siz,gca,0.9,0.7,0.7);
0172 graph1D_jd(0:dtn:tnmax,[0,normf0_1_1-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{2},2,siz,gca);
0173 graph1D_jd(0:dtn:tnmax,[0,(normf0_1_2-1)*1e3],0,0,'','','',leg,xlim,ylim,'-','none',colors{3},2,siz,gca);
0174 %
0175 set(gca,'xtick',[0:0.2:1]*tnmax)
0176 set(gca,'ytick',[-8:4:8]*1e-5)
0177 %
0178 figure(2),clf
0179 %
0180 leg = {'Free Mode','Fixed point','Normalization'};
0181 ylim = [-4,4]*1e-3;
0182 %
0183 graph1D_jd(0:dtn:tnmax,[0,normf0_2_0-1],0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none',colors{1},2,siz,gca,0.9,0.7,0.7);
0184 graph1D_jd(0:dtn:tnmax,[0,normf0_2_1-1],0,0,'','','',NaN,xlim,ylim,'-','none',colors{2},2,siz,gca);
0185 graph1D_jd(0:dtn:tnmax,[0,normf0_2_2-1],0,0,'','','',leg,xlim,ylim,'-','none',colors{3},2,siz,gca);
0186 %
0187 set(gca,'xtick',[0:0.2:1]*tnmax)
0188 set(gca,'ytick',[-4:2:4]*1e-3)
0189 %
0190 % residue
0191 %
0192 figure(3),clf
0193 %
0194 ylim = 10.^[-23,-21];
0195 ytick = 10.^[-23:-21];
0196 ylab = 'convergence parameter';
0197 %
0198 graph1D_jd(dtn:dtn:tnmax,residue_1_0,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','+',colors{1},2,siz,gca,0.9,0.7,0.7);
0199 graph1D_jd(dtn:dtn:tnmax,residue_1_1,0,1,'','','',NaN,xlim,ylim,'-','+',colors{2},2,siz,gca);
0200 graph1D_jd(dtn:dtn:tnmax,residue_1_2,0,1,'','','',leg,xlim,ylim,'-','+',colors{3},2,siz,gca);
0201 %
0202 set(gca,'xtick',[0:0.2:1]*tnmax)
0203 set(gca,'ytick',ytick)
0204 %
0205 figure(4),clf
0206 %
0207 graph1D_jd(dtn:dtn:tnmax,residue_2_0,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','+',colors{1},2,siz,gca,0.9,0.7,0.7);
0208 graph1D_jd(dtn:dtn:tnmax,residue_2_1,0,1,'','','',NaN,xlim,ylim,'-','+',colors{2},2,siz,gca);
0209 graph1D_jd(dtn:dtn:tnmax,residue_2_2,0,1,'','','',leg,xlim,ylim,'-','+',colors{3},2,siz,gca);
0210 %
0211 set(gca,'xtick',[0:0.2:1]*tnmax)
0212 set(gca,'ytick',ytick)
0213 %
0214 print_jd(2,'fig_thermal_norm_TiTe','./figures',1)
0215 print_jd(2,'fig_thermal_norm_TiTe2','./figures',2)
0216 %
0217 print_jd(2,'fig_thermal_norm_TiTe_res','./figures',3)
0218 print_jd(2,'fig_thermal_norm_TiTe2_res','./figures',4)
0219 %
0220 %************************************************************************************************************************************
0221 %
0222 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0223 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);