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 ***************************** 0019 % 0020 [qe,me,~,~,e0,~,~,mc2] = pc_dke_yp;%Universal physics constants 0021 % 0022 id_simul = 'Quench_ohm_nosync';%Simulation ID 0023 % 0024 % proc options 0025 % 0026 p_opt = 3;% printing and saving options : (-1) do nothing (0) print figures (1) print figures and save figures and results (2) save figures and results 0027 s_opt = true;% save distribution in file 0028 % 0029 f_opt = false;%true;% display electron distribution 0030 c_opt = true;% option for cylindrical calculation 0031 % 0032 % current and electric field 0033 % 0034 j0 = 1;%NaN;% initial current density (MA/m²) 0035 id_wave = '';%'LH_karney_mod_rev';%'';% for non empty id_wave, the initial current density is carried by the wave, with the diffusion coefficient adapted accordingly 0036 % 0037 opt.tol = 1e-3;% tolerance of electric field convergence 0038 opt.maxiconv = 20;% max # of iterations of electric field convergence 0039 opt.efield0 = 0.01;% initial guess for the electric field 0040 % 0041 % temperature evolution : 0042 % 0043 % o (0) : constant 0044 % o (1) : power law : Te = Te0*(1-t/tq)^(2/3) 0045 % o (2) : exp. law : Te = Tef + (Te0-Tef)*exp(-t/tq) 0046 % o (3) : step function : Te = Tef 0047 % o (4) : arena : Te = Tef + (Te0-Tef)*exp(-(t/tq)^2) 0048 % 0049 opt.T = 4;%1;% 0050 % 0051 betath0 = 0.1;% 0.01 corresponds to 5.1 keV 0052 betathf = 0.01;% 0.01 corresponds to 51 eV 0053 % 0054 % the ratio of the syncrotron reaction force to the collisional drag scales as wc2/wp2 0055 % -> set wpr = 0 or opt.synchro_mode = 0 to ignore syncrotron reaction force 0056 % 0057 wpr = 0;%([4]*qe/me)^2/([1e19]*qe^2/(me*e0));%wc2/wp2 0058 % 0059 opt.synchro_mode = 0; 0060 % 0061 % ion charge (single species, cold ions) 0062 % 0063 Zi = 1; 0064 % 0065 Emin = NaN;%runaway energy threshold in MeV (limit where electrons are "counted") 0066 % Emax = (sqrt(20^2*betath2+1)-1)*mc2/1e3;%1;%grid energy threshold in MeV [(sqrt(pnmax_S^2*betath2+1)-1)*mc2/1e3 = 0.63163 corresponds to pnmax_S = 20 for betath = 0.1] 0067 % 0068 % collision mode : 0069 % (0) : Relativistic Maxwellian background 0070 % (1) : High-velocity limit 0071 % (2) : Linearized Belaiev-Budker (momentum-conserving) 0072 % 0073 opt.coll_mode = 0; 0074 opt.bounce_mode = 0;% bounce-averaged calculation 0075 opt.boundary_mode_f = 0;% Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points 0076 opt.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 0077 % 0078 % time grid 0079 % 0080 % opt.tnq = 4;% opt.T = 0 or opt.T = 1 or opt.T = 3 0081 % opt.tnq = 2;% opt.T = 2 0082 opt.tnq = 50;% opt.T = 4 0083 % 0084 opt.tnmax = 500;%NaN;%5;%50;% 0085 % 0086 % if opt.tnmax = NaN, the simulation time tnmax is calculated as a function of opt.T and opt.tnq : 0087 % 0088 % - opt.T = 0 : opt.tnmax = opt.tnq 0089 % - opt.T = 1 : opt.tnmax = opt.tnq*(1 - (betathf/betath0)^3)) 0090 % - opt.T = 2 : opt.tnmax = 2*opt.tnq 0091 % - opt.T = 3 : opt.tnmax = opt.tnq 0092 % - opt.T = 4 : opt.tnmax = 2*opt.tnq 0093 % 0094 opt.nit = 100;%50;% number of time steps (lin scale) 0095 % 0096 graph.itdisp = 10:10:50;%[1,5,10];% display times - only used if opt.timevol > 0 0097 % 0098 % momentum grid 0099 % 0100 % momentum space grid -> you can specify nmhu_S and either 0101 % - the momentum grid pn_S directly 0102 % - the parameters np_S, pnmax_S, np_tail, pnmax_S_tail as follows 0103 % 0104 grid.nmhu_S = 101;% number of pitch angle grid points 0105 grid.np_S = 141;% number of bulk momentum grid points 0106 grid.pnmax_S = 28;%sqrt(((1 + Emax*1e3/mc2)^2 - 1)/betath2);% maximum of bulk momentum (normalized to pT) [pn = 28 corresponds to 1 MeV electrons for betath = 0.1] 0107 grid.np_tail = 331;%100;% number of tail momentum grid points 0108 grid.pnmax_S_tail = 597;% maximum of tail momentum (normalized to pT) [107,205,303,401,499,597] corresponds to [5,10,15,20,25,30] MeV electrons for betath = 0.1 0109 % 0110 % avalanche modelling : 0111 % 0112 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator 0113 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated 0114 % - For pnmin0_KO = NaN, the limit is set to pnmax_S 0115 % 0116 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator 0117 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy' 0118 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term 0119 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO 0120 % 0121 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator 0122 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used 0123 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference) 0124 % 0125 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator 0126 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons. 0127 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO 0128 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO 0129 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number" 0130 % 0131 opt.avalanche_mode = false;%true;% 0132 opt.pnmin0_KO = 28;% pnmin0_KO = 28 corresponds to 1 MeV electrons for betath = 0.1 0133 opt.pnmax1_KO = 0;%NaN;%grid.pnmax_S;% 0134 opt.pnmin2_KO = 1i;%NaN;%1;% 1i is the "natural" limit as for lower energies it is simply the identity operator 0135 opt.pnmax2_KO = NaN;%100;% 0136 % 0137 path_simul = '';%Simulation path 0138 % 0139 if opt.bounce_mode == 0, 0140 id_equil = 'TScyl';% cylindrical equilibrium, rho_S has no effect 0141 else 0142 id_equil = 'TScirc_e1';% toroidal equilibirum with a/R = 1 : then rho_S = r/R 0143 end 0144 % 0145 rho_S = 0.5; 0146 % 0147 opt.loadstructs = false;%true;% (true) use loadstructs, (false) recalculate all structures (for verification) 0148 % 0149 % graph.ylim_norm = [1e-35,1e1]; 0150 % graph.ytick_norm = 10.^(-35:5:5); 0151 % graph.ylim_gamma = [1e-34,1e-18]; 0152 % graph.ytick_gamma = 10.^(-34:4:-18); 0153 % graph.ylim_f = [1e-45,1e5]; 0154 % graph.ytick_f = 10.^(-45:5:5); 0155 % 0156 %************************************************************************************************************************************ 0157 % 0158 opt.pnmin_S = sqrt(((1 + Emin*1e3/mc2)^2 - 1))/betath0; 0159 % 0160 % RUN LUKE 0161 % 0162 [xTe,Efield,J,pfM,pf0,pfbM,savedata,trbetath] = frundke_quench(id_simul,j0,betath0,betathf,wpr,Zi,opt,grid,id_equil,id_wave,rho_S); 0163 % 0164 %************************************************************************************************************************************ 0165 % 0166 % PROCESS DATA 0167 % 0168 fproc_quench(xTe,Efield,J,pfM,pf0,pfbM,savedata,trbetath,graph,p_opt,f_opt,c_opt,id_simul); 0169 % 0170 %************************************************************************************************************************************ 0171 % 0172 if s_opt, 0173 savestr = [path_simul,'LUKEDIST_',id_simul,'.mat']; 0174 save(savestr,'-struct','savedata'); 0175 info_dke_yp(2,['Distribution data saved in ',savestr]); 0176 end 0177 % 0178 if p_opt > 0, 0179 filename = [path_simul,'DKE_RESULTS_',id_simul,'.mat']; 0180 save(filename); 0181 info_dke_yp(2,['Simulation data saved in ',filename]); 0182 end 0183 %