fprep_runaway

0001 function [dkepath,equil,dkeparam,dkedisplay,waves] = fprep_runaway(wpr,Zi,opt,grid,id_equil,id_wave,rho_S)
0002 %
0003 if nargin < 7,
0004    rho_S = 0.5;
0005 end
0006 %
0007 if nargin < 6,
0008    id_wave = '';
0009 end
0010 %
0011 if nargin < 4,
0012    grid = struct;
0013 end
0014 %
0015 if nargin < 3,
0016    opt = struct;
0017 end
0018 %
0019 if nargin < 2,
0020    Zi = 1;
0021 end
0022 %
0023 if nargin < 1,
0024    wpr = 0;
0025 end
0026 %
0027 if ~isfield(opt,'coll_mode'),
0028     opt.coll_mode = 0;
0029 end
0030 if ~isfield(opt,'bounce_mode'),
0031     opt.bounce_mode = 0;
0032 end
0033 if ~isfield(opt,'boundary_mode_f'),
0034     opt.boundary_mode_f = 0;
0035 end
0036 if ~isfield(opt,'norm_mode_f'),
0037     opt.norm_mode_f = 0;
0038 end
0039 if ~isfield(opt,'tn'),
0040     opt.tn = 10000;% typical time for asymptotic solution
0041 end
0042 if ~isfield(opt,'dtn'),
0043     opt.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn
0044 end
0045 if ~isfield(opt,'tnmin'),
0046     opt.tnmin = 0;
0047 end
0048 if ~isfield(opt,'synchro_mode'),
0049     if wpr > 0,
0050         opt.synchro_mode = 1;
0051     else
0052         opt.synchro_mode = 0;
0053     end
0054 end
0055 %
0056 if ~isfield(opt,'avalanche_mode'),
0057     opt.avalanche_mode = 0;
0058 end
0059 if ~isfield(opt,'pnmin0_KO'),% low threshold on initial primary electrons for the knock-on operator
0060     opt.pnmin0_KO = NaN;% NaN to enforce to pnmax_S
0061 end
0062 if ~isfield(opt,'pnmax1_KO'),% high threshold on initial secondary electrons for the knock-on operator
0063     opt.pnmax1_KO = 0;% NaN to enforce to pnmin2_KO
0064 end
0065 if ~isfield(opt,'pnmin2_KO'),% low threshold on final secondary electrons for the knock-on operator
0066     opt.pnmin2_KO = 1i;% 1i is the "natural" limit as for lower energies it is then merges with FP collisions
0067 end
0068 if ~isfield(opt,'pnmax2_KO'),% high threshold on final secondary electrons for the knock-on operator
0069     opt.pnmax2_KO = NaN;% NaN to enforce to pnmin0_KO
0070 end
0071 %
0072 if nargin < 5 || isempty(id_equil),
0073     if opt.bounce_mode == 0,
0074         id_equil = 'TScyl';%For plasma equilibrium
0075     else
0076         id_equil = 'TScirc_e1';%For plasma equilibrium
0077     end 
0078 end
0079 %
0080 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0081 %
0082 id_path = '';%For all paths used by DKE solver
0083 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0084 %
0085 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0086 %
0087 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0088 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0089 %
0090 id_display = 'NO_DISPLAY';%For output code display
0091 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0092 %
0093 path_wave = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0094 %
0095 %************************************************************************************************************************************
0096 %************************************************************************************************************************************
0097 %************************************************************************************************************************************
0098 %
0099 [dkepath,equil,dkeparam,dkedisplay,waves] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'waves',{id_wave},path_wave);
0100 %
0101 equil.pne = equil.pne/2;
0102 %
0103 %**************************************************************************
0104 %**************************************************************************
0105 %
0106 % luke parameters
0107 %
0108 if exist('dmumpsmex','file') == 3,
0109     dkeparam.invproc = -2;
0110 end
0111 %
0112 dkeparam.coll_mode = opt.coll_mode;% mode for the collision operator
0113 dkeparam.bounce_mode = opt.bounce_mode;% bounce averaging option
0114 %
0115 dkeparam.boundary_mode_f = opt.boundary_mode_f;%Enforcing the Maxwellian initial value at the first "boundary_mode_f" grid points
0116 dkeparam.norm_mode_f = opt.norm_mode_f;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0117 dkeparam.psin_S = psin_S;
0118 dkeparam.rho_S = rho_S;
0119 %
0120 if isfield(opt,'pnmin_S'),
0121     dkeparam.pnmin_S = opt.pnmin_S;
0122 end
0123 %
0124 % time grid
0125 %
0126 if opt.tnmin == 0,%linear time grid
0127     dkeparam.tn = opt.tn;% time grid (if array) or maximum time
0128     dkeparam.dtn = opt.dtn;% time grid (if array) or time step
0129 end
0130 if opt.tnmin > 0,%log time grid
0131     dkeparam.tn = logspace(log10(opt.tnmin),log10(opt.tn(1)),opt.dtn(1));
0132     if length(opt.tn) > 1,%log + lin
0133         tnlin = linspace(opt.tn(1),opt.tn(2),1+opt.dtn(2));        
0134         dkeparam.tn = [dkeparam.tn,tnlin(2:end)];
0135     end
0136     dkeparam.dtn = NaN;
0137 end
0138 %
0139 % momentum space grid
0140 %
0141 dkeparam.grid_uniformity = 0;
0142 if isfield(grid,'nmhu_S'),
0143     dkeparam.nmhu_S = grid.nmhu_S;
0144 end
0145 if isfield(grid,'sfac1'),
0146     dkeparam.sfac1 = grid.sfac1;
0147 end
0148 if isfield(grid,'snmhu1'),
0149     dkeparam.snmhu1 = grid.snmhu1;
0150 end
0151 if isfield(grid,'pn_S'), 
0152     dkeparam.pn_S = grid.pn_S;
0153     dkeparam.pnmax_S = max(grid.pn_S);%used for dNparmin, etc
0154 else
0155     if ~isfield(grid,'pnmax_S'),
0156         grid.pnmax_S = dkeparam.pnmax_S;
0157     end
0158     if ~isfield(grid,'np_S'),
0159         grid.np_S = dkeparam.np_S;
0160     end
0161     %
0162     pn_S = linspace(0,grid.pnmax_S,grid.np_S);
0163     %
0164     if ~isfield(grid,'np_bulk'),
0165         grid.np_bulk = 0;
0166     end
0167     %
0168     if grid.np_bulk == 0,
0169         pn_S_bulk = [];
0170     else
0171         pn_S = pn_S(pn_S > 1);
0172         pn_S_bulk = linspace(0,1,grid.np_bulk + 1);
0173     end
0174     %
0175     if ~isfield(grid,'np_tail') || ~isfield(grid,'pnmax_S_tail') || grid.pnmax_S_tail == 0,
0176         grid.np_tail = 0;
0177     end
0178     %
0179     if grid.np_tail == 0,
0180         pn_S_tail = [];
0181     else
0182         pn_S = pn_S(1:end-1);
0183         if imag(grid.pnmax_S_tail) > 0,
0184             grid.pnmax_S_tail = imag(grid.pnmax_S_tail)*grid.pnmax_S;
0185         end
0186         if isfield(grid,'opt') && grid.opt > 0,
0187             pn_S_tail = (linspace(grid.pnmax_S^(1/grid.opt),grid.pnmax_S_tail^(1/grid.opt),grid.np_tail + 1)).^grid.opt;
0188         else
0189             pn_S_tail = linspace(grid.pnmax_S,grid.pnmax_S_tail,grid.np_tail + 1);
0190         end
0191     end
0192     dkeparam.pn_S = [pn_S_bulk,pn_S,pn_S_tail];
0193     dkeparam.pnmax_S = grid.pnmax_S;%used for dNparmin, etc
0194 end
0195 %
0196 % avalanche modelling :
0197 %
0198 % dkeparam.pnmin0_KO is the low threshold on initial primary electrons for the knock-on operator
0199 % The condition pnmin0_KO <= pnmax_S is enforced so that the source can be calculated
0200 % - For pnmin0_KO = NaN, the limit is set to pnmax_S
0201 %
0202 % dkeparam.pnmax1_KO is the high threshold on initial secondary (bulk) electrons for the knock-on operator
0203 % in principle, pnmax1_KO <= pnmin2_KO so that secondary electrons can only 'lose energy'
0204 % - For pnmax1_KO = 0, the bulk population is taken as f(1)/fM(1), a logical estimate considering the sink term
0205 % - For pnmax1_KO = NaN, the limit is set to pnmin2_KO
0206 %
0207 % dkeparam.pnmin2_KO is the low threshold on final secondary electrons for the knock-on operator
0208 % - If pnmin2_KO == NaN, the current value of pc (or pnmin_S) is used
0209 % - If pnmin2_KO is imaginary, it is normalized to the actual thermal momentum (instead of reference)
0210 %
0211 % dkeparam.pnmax2_KO is the high threshold on final secondary electrons for the knock-on operator
0212 % in principle, as rosenbluth assumes p0 > p1 we must separate the domain of source runaways (ex : 1MeV) from that of secondary electrons.
0213 % Therefore, we should ensure pnmax2_KO <= pnmin0_KO
0214 % - For pnmax2_KO = NaN, the limit is set to pnmin0_KO
0215 % - If pnmax2_KO > pnmax_S (which violate the principle described above), the source between pnmax_S and pnmax2_KO is calculated as a "number"
0216 %
0217 dkeparam.avalanche_mode = opt.avalanche_mode;% knock on operator option
0218 dkeparam.pnmin0_KO = opt.pnmin0_KO;
0219 dkeparam.pnmax1_KO = opt.pnmax1_KO;
0220 dkeparam.pnmin2_KO = opt.pnmin2_KO;
0221 dkeparam.pnmax2_KO = opt.pnmax2_KO;
0222 %
0223 [qe,me,~,~,e0] = pc_dke_yp;%Universal physics constants
0224 %
0225 % synchrotron reaction force
0226 %
0227 dkeparam.synchro_mode = opt.synchro_mode;% synchrotron reaction force option
0228 if wpr > 0,
0229     if isfield(opt,'adjustB') && opt.adjustB, % adjust magnetic field amplitude  to have correct wp2/wc2 - synchrotron reaction force only
0230         pB = sqrt(equil.ptBx(:,1).^2 + equil.ptBy(:,1).^2 + equil.ptBPHI(:,1).^2).';
0231         pBfac = sqrt(equil.pne*qe^2*wpr/(me*e0))*me/qe./pB;
0232         ptBfac = repmat(pBfac.',[1,length(equil.theta)]);
0233         equil.ptBx = equil.ptBx.*ptBfac;
0234         equil.ptBy = equil.ptBy.*ptBfac;
0235         equil.ptBPHI = equil.ptBPHI.*ptBfac;
0236     else% adjust density to have correct wp2/wc2 - synchrotron reaction force only
0237         pB = sqrt(equil.ptBx(:,1).^2 + equil.ptBy(:,1).^2 + equil.ptBPHI(:,1).^2).';
0238         equil.pne = (pB*qe/me).^2*me*e0./(qe^2*wpr);
0239     end
0240 end
0241 %
0242 % equilibrium
0243 %
0244 equil.pzTi = 1e-10*ones(size(equil.pzTi));% cold ions
0245 %
0246 equil.zZi = [1,1,1,Zi];
0247 equil.pzni = [0;0;0;1/Zi]*equil.pne;% single ion species
0248 %
0249 if isfield(opt,'heavyions') && opt.heavyions,
0250     equil.zmi = 1e10*[1,1,1,1];
0251 end