rundke

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';%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 id_path = '';%For all paths used by DKE solver
0026 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0027 %
0028 psin_S = [0.001];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0029 rho_S = [];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0030 %
0031 id_equil = 'TScirc';%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 = '';%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 nit = 20;%number of iterations
0057 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)
0058 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
0059 %
0060 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0061 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0062 [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);
0063 dkeparam.boundary_mode_f = boundary_mode_f;
0064 dkeparam.coll_mode = coll_mode;
0065 dkeparam.dtn = dkeparam.tn;
0066 %
0067 dkeparam.psin_S = psin_S;
0068 dkeparam.rho_S = rho_S;
0069 %
0070 if exist('dmumpsmex');dkeparam.invproc = -2;end
0071 %
0072 [Znorm_u,Zcurr_u,ZP0_u,dke_out_u,radialDKE_u,equilDKE_u,momentumDKE_u,gridDKE_u,Zmomcoef_u,Zbouncecoef_u,Zmripple_u,mksa,XXsinksource_u] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0073 ne_u(1) = Znorm_u.x_0_fsav;
0074 for i = 2:nit,
0075    [Znorm_u,Zcurr_u,ZP0_u,dke_out_u] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out_u.XXf0,[]);
0076     ne_u(i) = Znorm_u.x_0_fsav;
0077 end
0078 %
0079 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0080 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0081 [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);
0082 dkeparam.boundary_mode_f = boundary_mode_f;
0083 dkeparam.coll_mode = coll_mode;
0084 dkeparam.dtn = dkeparam.tn;
0085 %
0086 dkeparam.psin_S = psin_S;
0087 dkeparam.rho_S = rho_S;
0088 %
0089 if exist('dmumpsmex');dkeparam.invproc = -2;end
0090 %
0091 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0092 ne(1) = Znorm.x_0_fsav;
0093 for i = 2:nit,
0094    [Znorm,Zcurr,ZP0,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out.XXf0,[]);
0095     ne(i) = Znorm.x_0_fsav;
0096 end
0097 %
0098 figure(1),clf,
0099 graph1D_jd([1:nit]*dkeparam.dtn,ne_u*mksa.ne_ref,0,0,'Integration time step (t_{coll})','Electron density (10^{+19} m^{-3)}',['Maxwellian distr., dtn = ',int2str(dkeparam.dtn)],NaN,'','','-','+','r',2,20,gca,1,0.6,0.6);hold on
0100 graph1D_jd([1:nit]*dkeparam.dtn,ne*mksa.ne_ref,0,0,'Integration time step (\tau_{coll})','Electron density (10^{+19} m^{-3)}',['Maxwellian distr., dtn = ',int2str(dkeparam.dtn)],NaN,'','','-','o','b',2,20,gca,1,0.6,0.6);
0101 axis([0,nit*dkeparam.dtn,mksa.ne_ref*0.98,mksa.ne_ref*1.02])
0102 line([0,20]*dkeparam.dtn,[mksa.ne_ref,mksa.ne_ref],'Color','red')
0103 legend('Uniform','Non-uniform','Location','SouthWest');
0104 %
0105 % Case with Ti = Te/2
0106 %
0107 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0108 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0109 [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);
0110 equil.pzTi = equil.pzTi/2;
0111 dkeparam.boundary_mode_f = boundary_mode_f;
0112 dkeparam.coll_mode = coll_mode;
0113 dkeparam.dtn = dkeparam.tn;
0114 %
0115 dkeparam.psin_S = psin_S;
0116 dkeparam.rho_S = rho_S;
0117 %
0118 if exist('dmumpsmex');dkeparam.invproc = -2;end
0119 %
0120 [Znorm_u1,Zcurr_u1,ZP0_u1,dke_out_u1,radialDKE_u1,equilDKE_u1,momentumDKE_u1,gridDKE_u1,Zmomcoef_u1,Zbouncecoef_u1,Zmripple_u1,mksa,XXsinksource_u1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0121 ne_u1(1) = Znorm_u1.x_0_fsav;
0122 for i = 2:nit,
0123    [Znorm_u1,Zcurr_u1,ZP0_u1,dke_out_u1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out_u1.XXf0,[]);
0124     ne_u1(i) = Znorm_u1.x_0_fsav;
0125 end
0126 %
0127 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0128 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0129 [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);
0130 equil.pzTi = equil.pzTi/2;
0131 dkeparam.boundary_mode_f = boundary_mode_f;
0132 dkeparam.coll_mode = coll_mode;
0133 dkeparam.dtn = dkeparam.tn;
0134 %
0135 dkeparam.psin_S = psin_S;
0136 dkeparam.rho_S = rho_S;
0137 %
0138 if exist('dmumpsmex');dkeparam.invproc = -2;end
0139 %
0140 [Znorm1,Zcurr1,ZP01,dke_out1,radialDKE1,equilDKE1,momentumDKE1,gridDKE1,Zmomcoef1,Zbouncecoef1,Zmripple1,mksa,XXsinksource1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0141 ne1(1) = Znorm1.x_0_fsav;
0142 for i = 2:nit,
0143    [Znorm1,Zcurr1,ZP01,dke_out1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out1.XXf0,[]);
0144     ne1(i) = Znorm1.x_0_fsav;
0145 end
0146 %
0147 figure(2),clf,graph1D_jd([1:nit]*dkeparam.dtn,ne_u1*mksa.ne_ref,0,0,'Integration time step (t_{coll})','Electron density (10^{+19} m^{-3)}',['Maxwellian distr., dtn = ',int2str(dkeparam.dtn)],NaN,'','','-','+','r',2,20,gca,1,0.6,0.6);
0148 graph1D_jd([1:nit]*dkeparam.dtn,ne1*mksa.ne_ref,0,0,'Integration time step (\tau_{coll})','Electron density (10^{+19} m^{-3)}',['Ti = Te/2, Maxwellian distr., dtn = ',int2str(dkeparam.dtn)],NaN,'','','-','o','b',2,20,gca,1,0.6,0.6);
0149 axis([0,nit*dkeparam.dtn,mksa.ne_ref*0.97,mksa.ne_ref*1.03])
0150 line([0,20]*dkeparam.dtn,[mksa.ne_ref,mksa.ne_ref],'Color','red')
0151 legend('Uniform','Non-uniform','Location','SouthWest');
0152 %
0153 print_jd(p_opt,'thermaldensity','./figures',1)
0154 %
0155 %************************************************************************************************************************************
0156 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0157 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);