frundke_quench

0001 function [xTe,Efield,J,pfM,pf0,pfbM,savedata,trbetath] = frundke_quench(id_simul,j0,betath0,betathf,wpr,Zi,opt,varargin)
0002 %
0003 [dkepath,equil0,dkeparam,dkedisplay,waves] = fprep_runaway(wpr,Zi,opt,varargin{:});
0004 %
0005 [~,~,~,~,~,~,~,mc2] = pc_dke_yp;%Universal physics constants
0006 %
0007 %************************************************************************************************************************************
0008 %
0009 % initial (time-asymptotic) distribution
0010 %
0011 equil0.pTe = betath0^2*mc2*ones(size(equil0.pTe));%initial temperature
0012 %
0013 dkeparam0 = dkeparam;% parameters for initial distribution calculation
0014 %
0015 dkeparam0.tn = NaN;%time-asymptotic
0016 dkeparam0.dtn = NaN;%time-asymptotic
0017 dkeparam0.opt_struct = 2;%to keep all fields in output structures for loadstruct
0018 %
0019 if isempty(waves),% the plasma is initially in ohmic state
0020     %
0021     dkeparam0.Jtot_ohm = j0;% current density
0022     %
0023     dkeparam0.nit_ohm = opt.maxiconv;
0024     dkeparam0.tol_ohm = opt.tol;
0025     dkeparam0.efield_ohm = opt.efield0;
0026     %
0027     [Znorm0,Zcurr0,ZP00,dke_out0,radialDKE0,equilDKE0,momentumDKE0,gridDKE0,Zmomcoef0,Zbouncecoef0,Zmripple0,mksa0,XXsinksource0] = main_dke_yp(id_simul,dkepath,equil0,dkeparam0,dkedisplay,'','','','',[],[]);
0028     %
0029 else% the plasma is initially in RF driven state
0030     %
0031     for it_wave = 1:opt.maxiconv,
0032         %
0033         [Znorm0,Zcurr0,ZP00,dke_out0,radialDKE0,equilDKE0,momentumDKE0,gridDKE0,Zmomcoef0,Zbouncecoef0,Zmripple0,mksa0,XXsinksource0] = main_dke_yp(id_simul,dkepath,equil0,dkeparam0,dkedisplay,'',waves,'','',[],[]);
0034         %
0035         Jtot_rf = Zcurr0.x_0_fsav*mksa0.j_ref;
0036         %
0037         if abs(Jtot_rf - j0)/j0 <= opt.tol,
0038             disp(['convergence obtained for RF current after ',num2str(it_wave),' iterations']);
0039             break
0040         elseif it_wave == opt.maxiconv,
0041             disp(['WARNING : convergence not obtained for RF current after ',num2str(it_wave),' iterations']);
0042         else
0043             waves{1}.rays{1}.P0_2piRp = waves{1}.rays{1}.P0_2piRp*j0/Jtot_rf;
0044         end
0045         %
0046     end
0047     %
0048 end
0049 %
0050 dkeparam.Jtot_ohm = j0;% current density
0051 %
0052 dkeparam.nit_ohm = opt.maxiconv;
0053 dkeparam.tol_ohm = opt.tol;
0054 dkeparam.efield_ohm = opt.efield0;
0055 %
0056 dkeparam.opt_struct = 0;%to keep all fields in output structures for loadstruct
0057 %
0058 dmhu = momentumDKE0.dmhu;
0059 mhu = momentumDKE0.mhu;
0060 dpn = momentumDKE0.dpn;
0061 pn = momentumDKE0.pn;
0062 pn2 = gridDKE0.pn2;
0063 gamma = Zmomcoef0.gamma;
0064 %
0065 if opt.loadstructs,
0066     loadstructs.radialgrid_dke_jd.radialDKE = radialDKE0;
0067     %
0068 %     loadstructs.equilibrium_jd.equilDKE = equilDKE0;
0069     %
0070     loadstructs.momentumgrid_dke_jd.momentumDKE = momentumDKE0;
0071     loadstructs.matrix3Dgridbuilder_dke_yp.gridDKE = gridDKE0;
0072     loadstructs.momentumcoefbuilder_dke_yp.Zmomcoef = Zmomcoef0;
0073     loadstructs.bouncecoefbuilder_dke_yp.Zbouncecoef = Zbouncecoef0;
0074     loadstructs.mripplecoefbuilder_dke_yp.Zmripple = Zmripple0;
0075     loadstructs.mksacoefbuilder_dke_yp.mksa = mksa0;
0076     loadstructs.sinksourcecoeffbuilder_dke_yp.XXsinksource = XXsinksource0;
0077     %
0078     loadstructs.rad_dke_yp.ZXXD_a = dke_out0.ZXXD_a;
0079     loadstructs.rad_dke_yp.ZXXF_a = dke_out0.ZXXF_a;
0080     %
0081 %     loadstructs.coll_dke_jd.ZXXD_c = dke_out0.ZXXD_c;
0082 %     loadstructs.coll_dke_jd.ZXXF_c = dke_out0.ZXXF_c;
0083 %     loadstructs.coll_dke_jd.ZXXD_c_tp = NaN;%dke_out.ZXXD_c_tp;
0084 %     loadstructs.coll_dke_jd.ZXXF_c_tp = NaN;%dke_out.ZXXF_c_tp;
0085 %     loadstructs.coll_dke_jd.XXfM = dke_out0.XXfM;
0086 %     loadstructs.coll_dke_jd.XXILor = dke_out0.XXILor;
0087     %
0088     loadstructs.efield_dke_jd.ZXXD_e = dke_out0.ZXXD_e;
0089     loadstructs.efield_dke_jd.ZXXF_e = dke_out0.ZXXF_e;
0090     loadstructs.efield_dke_jd.ZXXD_e_tp = NaN;%dke_out.ZXXD_e_tp;
0091     loadstructs.efield_dke_jd.ZXXF_e_tp = NaN;%dke_out.ZXXF_e_tp;
0092     loadstructs.efield_dke_jd.xepsi_init = dke_out0.xepsi_init;
0093     loadstructs.efield_dke_jd.xEfield_validity = dke_out0.xEfield_validity;
0094     %
0095 %     loadstructs.rfwave_dke_jd.waveparam = dke_out0.waveparam;
0096 %     loadstructs.rfwave_dke_jd.xyprop_dke = dke_out0.xyprop_dke;
0097 %     loadstructs.rfwave_dke_jd.xyP0_2piRp_mod = dke_out0.xyP0_2piRp_mod;
0098 %     loadstructs.rfwave_dke_jd.yb = dke_out0.yb;
0099 %     loadstructs.rfwave_dke_jd.yP0_2piRp = dke_out0.yP0_2piRp;
0100 %     loadstructs.rfwave_dke_jd.xys = dke_out0.xys;
0101     %
0102 %     loadstructs.loop_rfdiff_dke_jd.ZXYD_rf = dke_out0.ZXYD_rf;
0103 %     loadstructs.loop_rfdiff_dke_jd.ZXYF_rf = dke_out0.ZXYF_rf;
0104 %     loadstructs.loop_rfdiff_dke_jd.ZXYD_rf_tp = NaN;%dke_out.ZXYD_rf_tp;
0105 %     loadstructs.loop_rfdiff_dke_jd.ZXYF_rf_tp = NaN;%dke_out.ZXYF_rf_tp;
0106 %     loadstructs.loop_rfdiff_dke_jd.gridindex_rf = dke_out0.gridindex_rf;
0107     %
0108     % loadstructs.finit_dke_yp.XXfinit = dke_out0.XXfinit;
0109     %
0110 else
0111     %
0112     loadstructs = [];
0113     %
0114     dkeparam.Te_ref = mksa0.Te_ref;%Reference temperature
0115     dkeparam.ne_ref = mksa0.ne_ref;%Reference density
0116     %
0117 end
0118 %
0119 I_tot = Zcurr0.I_tot;
0120 %
0121 XXfinit = dke_out0.XXf0;
0122 %
0123 ohm = struct;
0124 ohm.id = 'evolution';
0125 ohm.rho = gridDKE0.xrho;
0126 ohm.epsi = dke_out0.xepsi{end}*mksa0.Edreicer_ref;% guess is previous time
0127 ohm.xnrem_init = dke_out0.xnorm_rem{end};
0128 ohm.xnrep_init = dke_out0.xnorm_rep{end};
0129 %
0130 %************************************************************************************************************************************
0131 %
0132 % temperature evolution and time dependent distribution
0133 %
0134 % time grid
0135 %
0136 nit = opt.nit;
0137 %
0138 if isnan(opt.tnmax),
0139     if opt.T == 0,
0140         opt.tnmax = opt.tnq;
0141     elseif opt.T == 1,
0142         opt.tnmax = opt.tnq*(1 - (betathf/betath0)^3);
0143     elseif opt.T == 2,
0144         opt.tnmax = 2*opt.tnq;
0145     elseif opt.T == 3,
0146         opt.tnmax = opt.tnq;
0147     elseif opt.T == 4,
0148         opt.tnmax = 2*opt.tnq;
0149     end
0150 end
0151 %
0152 % time grid (normalized to collision time) -> you can specify :
0153 % - linear arbitrary grid with tn array                 : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0
0154 % - linear arbitrary grid with dtn array                : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10);
0155 % - linear grid with tnmax and grid step                : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000;
0156 % - linear grid with tnmax and number of grid steps     : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i;
0157 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10;
0158 %
0159 dkeparam.tn = NaN;
0160 dkeparam.dtn = opt.tnmax/nit;
0161 tn = linspace(dkeparam.dtn,opt.tnmax,nit);
0162 %
0163 if opt.T == 0,
0164     tbetath = repmat(betath0,[1,nit]);
0165 elseif opt.T == 1,
0166     tbetath = betath0*(1-tn/opt.tnq).^(1/3); 
0167 elseif opt.T == 2,
0168     tbetath = sqrt(betathf^2 + (betath0^2-betathf^2)*exp(-tn/opt.tnq)); 
0169 elseif opt.T == 3,
0170     tbetath = repmat(betathf,[1,nit]); 
0171 elseif opt.T == 4,
0172     tbetath = sqrt(betathf^2 + (betath0^2-betathf^2)*exp(-(tn/opt.tnq).^2)); 
0173 end
0174 %
0175 trbetath = [1,tbetath/betath0];
0176 %
0177 tequil = repmat(equil0,[1,nit]);
0178 %
0179 for it = 1:nit,
0180     tequil(it).pTe = tbetath(it)^2*mc2*ones(size(tequil(it).pTe));
0181 end
0182 %
0183 for it = 1:nit,    
0184     [Znorm(it),Zcurr(it),ZP0(it),dke_out(it)] = main_dke_yp(id_simul,dkepath,tequil(it),dkeparam,dkedisplay,ohm,'','','',XXfinit{it},[],loadstructs);
0185     XXfinit{it+1} = dke_out(it).XXf0{1};
0186     %
0187     ohm.epsi = dke_out(it).xepsi{end}*mksa0.Edreicer_ref;% guess is previous time
0188     ohm.xnrem_init = dke_out(it).xnorm_rem{end};
0189     ohm.xnrep_init = dke_out(it).xnorm_rep{end};
0190     %
0191     disp(['Time step ',num2str(it),'/',num2str(nit),' calculated.'])
0192 end
0193 %
0194 % calculation of the Spitzer current
0195 %
0196 dkeparam.avalanche_mode = 0;
0197 dkeparam.dtn = NaN;%time-asymptotic
0198 dkeparam.nit_ohm = 1;%no field adjustment
0199 %dkeparam = rmfield(dkeparam,{'Te_ref','ne_ref'});%TODO : pb when we keep the old grid...
0200 %
0201 ohm = struct;
0202 ohm.id = 'linear';
0203 ohm.rho = gridDKE0.xrho;
0204 ohm.epsi = 1e-6;% linear limit
0205 %
0206 if opt.loadstructs,
0207     loadstructs = rmfield(loadstructs,'efield_dke_jd');
0208 end
0209 %
0210 [~,Zcurr_lin0] = main_dke_yp(id_simul,dkepath,equil0,dkeparam,dkedisplay,ohm,'','','',[],[],loadstructs);
0211 for it = 1:nit,    
0212     [~,Zcurr_lin(it),~,~,~,~,~,~,~,~,~,mksa_lin(it)] = main_dke_yp(id_simul,dkepath,tequil(it),dkeparam,dkedisplay,ohm,'','','',[],[],loadstructs);
0213     disp(['Time step ',num2str(it),'/',num2str(nit),' calculated.'])
0214 end
0215 %
0216 pfM = cell(1,nit+1);
0217 pf0 = cell(1,nit+1);
0218 %
0219 xTe.normM = NaN(1,nit+1);
0220 xTe.norm0 = NaN(1,nit+1);
0221 %
0222 pfbM = cell(1,nit+1);
0223 xTe.normb = NaN(1,nit+1);
0224 xTe.norms = NaN(1,nit+1);
0225 xTe.delta = NaN(1,nit+1);
0226 err = NaN(1,nit+1);
0227 %
0228 pfM{1} = (1/2)*integral_dke_jd(dmhu,dke_out0.XXfM,2);%mhu-averaged distribution
0229 pf0{1} = (1/2)*integral_dke_jd(dmhu,dke_out0.XXf0{1},2);%mhu-averaged distribution
0230 %
0231 for it = 1:nit,
0232     pfM{it+1} = (1/2)*integral_dke_jd(dmhu,dke_out(it).XXfM,2);%mhu-averaged distribution
0233     pf0{it+1} = (1/2)*integral_dke_jd(dmhu,dke_out(it).XXf0{1},2);%mhu-averaged distribution
0234 end
0235 %
0236 for it = 1:nit+1,
0237     %
0238     xTe.normM(it) = 4*pi*integral_dke_jd(dpn.',pn2.'.*pn2.'.*pfM{it}./gamma.',1)/3;
0239     xTe.norm0(it) = 4*pi*integral_dke_jd(dpn.',pn2.'.*pn2.'.*pf0{it}./gamma.',1)/3;
0240     %
0241     [xTe.normb(it),xTe.norms(it),xTe.delta(it),pfbM{it},err(it)] = bimaxfit_jd(dpn.',pn.',gamma.',pf0{it},1,3,xTe.normM(it),1);
0242 end
0243 %
0244 J.tot = NaN(1,nit+1);% total current density
0245 J.lin = NaN(1,nit+1);% Spitzer current density
0246 J.rei = NaN(1,nit+1);% Internal RE current density
0247 J.ree = NaN(1,nit+1);% External RE current density
0248 %
0249 Efield = NaN(1,nit+1);%E/E_c
0250 %
0251 Efield(1) = dke_out0.xepsi{end}/mksa0.betath_ref^2;
0252 %
0253 J.tot(1) = dke_out0.xcurr{end};% total current density
0254 J.lin(1) = Zcurr_lin0.x_0_fsav*(dke_out0.xepsi{end}.*mksa0.Edreicer_ref)/ohm.epsi;% Spitzer current density
0255 J.rei(1) = Zcurr0.x_Rfrac_0_fsav;% Internal RE current density
0256 J.ree(1) = dke_out0.xcurr_rem{end} + dke_out0.xcurr_rep{end};% External RE current density
0257 %
0258 for it = 1:nit,
0259     Efield(it+1) = dke_out(it).xepsi{1}/mksa0.betath_ref^2;
0260     %
0261     J.tot(it+1) = dke_out(it).xcurr{end};% total current density
0262     %
0263     % maybe we must multiply Zcurr_lin(it) by the bulk density xnorm0 or ((1-delta)*xnorm0 ?
0264     J.lin(it+1) = Zcurr_lin(it).x_0_fsav*(mksa_lin(it).j_ref/mksa0.j_ref)*(dke_out(it).xepsi{end}*mksa0.Edreicer_ref)/ohm.epsi;% Spitzer current density
0265     J.rei(it+1) = Zcurr(it).x_Rfrac_0_fsav;% Internal RE current density
0266     J.ree(it+1) = dke_out(it).xcurr_rem{end} + dke_out(it).xcurr_rep{end};% External RE current density
0267 end
0268 %
0269 J.hot = J.tot - J.lin - J.rei - J.ree;% Hot tail current density
0270 %
0271 J.j0 = j0/mksa0.j_ref;   
0272 %
0273 [np,nmhu] = size(dke_out0.XXf0{1});
0274 %
0275 savedata = struct;
0276 %
0277 savedata.t = [0,tn];
0278 savedata.betath0 = betath0;
0279 savedata.pn = pn;
0280 savedata.dpn = dpn;
0281 savedata.xi = mhu;
0282 savedata.dxi = dmhu;
0283 %
0284 savedata.rbetath = [betath0,tbetath]/betath0;
0285 savedata.f = NaN(np,nmhu,nit+1);
0286 savedata.fM = NaN(np,nmhu,nit+1);
0287 savedata.f(:,:,1) = dke_out0.XXf0{1}/Znorm0.x_0_fsav;
0288 savedata.fM = dke_out0.XXfM;
0289 for it = 1:nit,
0290     savedata.f(:,:,it+1) = dke_out(it).XXf0{1}/Znorm(it).x_0_fsav;
0291     savedata.fM(:,:,it+1) = dke_out(it).XXfM;
0292 end
0293 %
0294 savedata.info.t = 'time grid, normalized to electron thermal collision time 4*pi*e0^2*m^2*vT0^3/(ne*e^4)';
0295 savedata.info.betath0 = 'ratio of the reference thermal velocity VT0=sqrt(T0/m) to the speed of light c. Also equal to sqrt(T0/mc2)';
0296 savedata.info.pn = 'momentum grid, normalized to the reference thermal momentum pT0 = m*vT0';
0297 savedata.info.dpn = 'momentum grid steps';
0298 savedata.info.xi = 'pitch angle grid, xi = p||/p';
0299 savedata.info.dxi = 'pitch angle grid steps';
0300 savedata.info.rbetath = 'temperature evolution prescribed through field particles in the collision operator: rbetath = vT/vT0';
0301 savedata.info.f = 'electron distribution function as a function of (pn,xi,t), normalized to the electron density';
0302 savedata.info.fM = 'Maxwellian distribution function as a function of (pn,xi,t), normalized to the electron density';
0303 %