Description of frundke

0001 function [outputs] = frundke_metis(simul,opt,grid,initstate,metisoptfile,metisoptions)
0002 %
0003 dkepath = load_structures_yp('dkepath','','');
0004 %
0005 % initialization of the metis structure
0006 %
0007 [option,cons1t,geo1t] = zerodevolution;
0008 %
0009 if nargin > 4 && ~isempty(metisoptfile),
0010     metisoptdata = load(metisoptfile);
0011     option = conc_struct_jd(option,metisoptdata.z0dinput.option);% update options with specific set from METIS file
0012 end
0013 %
0014 if nargin > 5 && isstruct(metisoptions),% update options with user input structure
0015     option = conc_struct_jd(option,metisoptions);
0016 end
0017 %
0018 % first time is set by initial distribution
0019 %
0020 equil = initstate.luke_input.equil;% LUKE equilibrium
0021 wavestructs = initstate.luke_input.wavestructs;% LUKE input wave structures
0022 %
0023 geo1t.a = (max(equil.ptx(end,:)) - min(equil.ptx(end,:)))/2;
0024 geo1t.R = equil.Rp;
0025 geo1t.K = equil.pelong(end);
0026 geo1t.d = equil.ptrian(end);
0027 geo1t.z0 = equil.Zp;
0028 geo1t.b0 = equil.ptBPHI(end,1)*(1 + equil.ptx(end,1)/equil.Rp);% based on constant magnetic rigidity
0029 geo1t.vp = [];% calculated by METIS
0030 geo1t.sp = [];% calculated by METIS
0031 geo1t.sext = [];% calculated by METIS
0032 %
0033 cons1t.temps = 0;
0034 cons1t.ip = equil.ip;
0035 cons1t.flux = 0;% ip used instead
0036 %
0037 pneh = (equil.pne(1:end-1) + equil.pne(2:end)).'/2;
0038 pdxh = diff(equil.ptx(:,1));
0039 cons1t.nbar = sum(pneh.*pdxh)/equil.ptx(end,1);
0040 %
0041 % line averaged density
0042 %
0043 pzeff = (equil.zZi.^2*equil.pzni)./equil.pne;
0044 pzeffh = (pzeff(1:end-1) + pzeff(2:end)).'/2;
0045 cons1t.zeff = sum(pzeffh.*pdxh)/equil.ptx(end,1);
0046 %
0047 % line averaged zeff
0048 %
0049 cons1t.picrh = 0;
0050 cons1t.pnbi = 0;
0051 cons1t.ftnbi = 0;
0052 %
0053 plh = 0;
0054 wpecrh = zeros(1,0);
0055 wwecrh = zeros(1,0);
0056 %
0057 if opt.waves,% initial current is associated with RF waves
0058     for wavestruct = wavestructs,% sum ofver wave input structures
0059         %
0060         launch = wavestruct{1}.launch;
0061         %
0062         if strcmp(launch.type,'LH'),
0063             plh = plh + sum(launch.bPlhtot);
0064         elseif strcmp(launch.type,'EC'), 
0065             wpecrh(1,end+1) = sum(launch.yP_L);
0066             wwecrh(1,end+1) = launch.omega_rf;
0067         else
0068             disp(['wave type ',launch.type,' not recognized. Discarded.'])
0069         end
0070     end
0071 end
0072 cons1t.plh = plh;
0073 cons1t.pecrh = sum(wpecrh);
0074 %
0075 % ECRH deposition
0076 %
0077 n = [1,2,3];% harmonic numbers
0078 [qe,me] = pc_dke_yp;%Universal physics constants
0079 %
0080 nwBres = (1./n.')*wwecrh*me/qe;% resonant B field
0081 nwRres = geo1t.R*geo1t.b0./nwBres;% resonant major radius (in vacuum)
0082 nwabsxres = abs(nwRres - geo1t.R);% resonant |R-Rp|
0083 wabsxres = min([ones(1,length(wwecrh));nwabsxres]);% minimum resonant |R-Rp| (versus harmonic number)
0084 if cons1t.pecrh > 0,    
0085     xece = sum(wabsxres.*wpecrh)/cons1t.pecrh;
0086 else
0087     xece = 0;
0088 end
0089 %
0090 cons1t.hmore = 1;% confinement factor adjustment
0091 cons1t.iso = 0;% T fraction in DT mix
0092 %
0093 sepa1t.Rsepa = equil.Rp + equil.ptx(end,:);
0094 sepa1t.Zsepa = equil.Zp + equil.pty(end,:);
0095 %
0096 % initialization of z0dstruct
0097 %
0098 if isappdata(0,'NE_EXP'),
0099     disp('Clearing external electron density data')
0100     rmappdata(0,'NE_EXP');
0101 end
0102 if isappdata(0,'TE_EXP'),
0103     disp('Clearing external electron temperature data')
0104     rmappdata(0,'TE_EXP');
0105 end
0106 if isappdata(0,'LH_SHAPE_EXP'),
0107     disp('Clearing external LH wave data')
0108     rmappdata(0,'LH_SHAPE_EXP');
0109 end
0110 if isappdata(0,'ECCD_SHAPE_EXP'),
0111     disp('Clearing external EC wave data')
0112     rmappdata(0,'ECCD_SHAPE_EXP');
0113 end
0114 if isappdata(0,'RUNAWAY_EXP'),
0115     disp('Clearing external runaway data')
0116     rmappdata(0,'RUNAWAY_EXP');
0117 end
0118 %
0119 [zs_init,profil_init,z0dstruct_init] = zerodevolution([],option,0,cons1t,geo1t,[],sepa1t,[]);
0120 %
0121 % external data
0122 %
0123 if opt.necons,
0124     NE_EXP.temps = [-1,0];
0125     NE_EXP.x = equil.pli;
0126     NE_EXP.ne = [1;1]*equil.pne;
0127     %
0128     disp('Loading external electron density data')
0129     setappdata(0,'NE_EXP',NE_EXP)
0130 end
0131 %
0132 if opt.tecons,
0133     TE_EXP.temps = [-1,0];
0134     TE_EXP.x = equil.pli;
0135     TE_EXP.te = [1;1]*equil.pTe*1e3;
0136     %
0137     disp('Loading external electron temperature data')
0138     setappdata(0,'TE_EXP',TE_EXP)
0139 end
0140 %
0141 if opt.jcons > 0,
0142     equilDKE = initstate.equilDKE;
0143     Zcurr = initstate.Zcurr;
0144     ZP0 = initstate.ZP0;
0145     mksa = initstate.mksa;
0146     %
0147     prho = equil.ptx(:,1)/equil.ptx(end,1);
0148     xli = interp1(prho,equil.pli,equilDKE.xrho,'linear');
0149     jluke = mksa.j_ref*1e6*Zcurr.x_0_vcfsav;
0150     prf = mksa.P_ref*1e6*ZP0.x_rf_fsav.';
0151     %
0152     if isstruct(initstate.luke_input.ohm),
0153         if isfield(initstate,'Zcurr_lin') && isfield(initstate,'ohm_lin'),% Spitzer contribution is calculated
0154             %
0155             Zcurr_lin = initstate.Zcurr_lin;
0156             dke_out_lin = initstate.Zcurr_lin;
0157             %
0158         else% Spitzer contribution not calculated
0159             %
0160             ohm_lin = initstate.luke_input.ohm;
0161             ohm_lin.epsi = 1e-6*ones(size(ohm_lin.epsi));% linear limit for Spitzer calculation
0162             %
0163             [~,Zcurr_lin,~,dke_out_lin] = main_dke_yp('lin',dkepath,initstate.luke_input.equil,initstate.luke_input.dkeparam,initstate.luke_input.dkedisplay,ohm_lin,'','','',[],[]);
0164             %
0165         end
0166         %
0167         dke_out = initstate.dke_out;
0168         %
0169         johm = mksa.j_ref*1e6*Zcurr_lin.x_0_vcfsav.*dke_out.xepsi./dke_out_lin.xepsi;
0170         %
0171     else
0172         johm = zeros(size(xli));
0173     end        
0174     %
0175     RUNAWAY_EXP.jrun = zeros(2,length(xli));
0176     LH_SHAPE_EXP.jlh = zeros(2,length(xli));
0177     ECCD_SHAPE_EXP.jeccd = zeros(2,length(xli));
0178     %
0179     if cons1t.plh > 0 && cons1t.pecrh > 0,
0180         disp(['current prescription not implemented with several wave types combined. Discarded.'])
0181     elseif cons1t.plh > 0
0182         %
0183         LH_SHAPE_EXP.temps = [-1,0];
0184         LH_SHAPE_EXP.x = xli;
0185         LH_SHAPE_EXP.jlh = [1;1]*(jluke - johm);
0186         LH_SHAPE_EXP.plh = [1;1]*prf;
0187         %
0188         if opt.jcons > 1,
0189             LH_SHAPE_EXP.nbar = [1;1]*cons1t.nbar;% to account for density variations in eta_lh
0190         end
0191         %
0192         disp('Loading external LH wave data')
0193         setappdata(0,'LH_SHAPE_EXP',LH_SHAPE_EXP)
0194         %
0195     elseif cons1t.pecrh > 0
0196         %
0197         ECCD_SHAPE_EXP.temps = [-1,0];
0198         ECCD_SHAPE_EXP.x = xli;
0199         ECCD_SHAPE_EXP.jeccd = [1;1]*(jluke - johm);
0200         ECCD_SHAPE_EXP.peccd = [1;1]*prf;
0201         %
0202         if opt.jcons > 1,
0203             ECCD_SHAPE_EXP.nbar = [1;1]*cons1t.nbar;% to account for density variations in eta_ec
0204         end
0205         %
0206         disp('Loading external EC wave data')
0207         setappdata(0,'ECCD_SHAPE_EXP',ECCD_SHAPE_EXP)
0208         %
0209     else% initial LUKE current was all runaway
0210         %
0211         RUNAWAY_EXP.temps = [-1,0];
0212         RUNAWAY_EXP.x = xli;
0213         RUNAWAY_EXP.jrun = [1;1]*(jluke - johm);
0214         RUNAWAY_EXP.irun = 1e6*[1;1]*Zcurr.I_tot;
0215         RUNAWAY_EXP.iohm = [1;1]*NaN;% for first order nonlinear corrections - not used yet
0216         RUNAWAY_EXP.ip = [1;1]*cons1t.ip;
0217         %
0218         disp('Loading external runaway data')
0219         setappdata(0,'RUNAWAY_EXP',RUNAWAY_EXP)
0220         %
0221     end
0222 end
0223 %
0224 % initial METIS state calculation with external data
0225 %
0226 [zs,profil,z0dstruct] = zerodevolution(z0dstruct_init,option,0,cons1t,geo1t,[],sepa1t,[]);
0227 %
0228 % adjust initial prf for non-inductive state
0229 %
0230 if opt.ni,
0231     if cons1t.plh > 0 && cons1t.pecrh > 0,
0232         disp(['power/current adjustment not implemented with several wave types combined. Discarded.'])
0233     elseif cons1t.plh > 0,
0234         %
0235         LH_SHAPE_EXP.jlh = LH_SHAPE_EXP.jlh*(1 + zs.iohm/zs.ilh);
0236         disp('Loading external LH wave data')
0237         setappdata(0,'LH_SHAPE_EXP',LH_SHAPE_EXP)
0238         %
0239 %         cons1t.plh = cons1t.plh*(1 + zs.iohm/zs.ilh);
0240         %
0241     elseif cons1t.pecrh > 0,
0242         %
0243         ECCD_SHAPE_EXP.jeccd = ECCD_SHAPE_EXP.jeccd*(1 + zs.iohm/zs.ieccd);
0244         disp('Loading external EC wave data')
0245         setappdata(0,'ECCD_SHAPE_EXP',ECCD_SHAPE_EXP)
0246         %
0247 %         cons1t.peccd = cons1t.peccd*(1 + zs.iohm/zs.ieccd);
0248         %
0249     else
0250         %
0251         RUNAWAY_EXP.jrun = RUNAWAY_EXP.jrun*(1 + zs.iohm/zs.irun);
0252         RUNAWAY_EXP.irun = [1;1]*zs.irun*(1 + zs.iohm/zs.irun);
0253         disp('Loading external runaway data')
0254         setappdata(0,'RUNAWAY_EXP',RUNAWAY_EXP)
0255         %
0256     end
0257     %
0258     [zs,profil,z0dstruct] = zerodevolution(z0dstruct,option,0,cons1t,geo1t,[],sepa1t,[]);
0259     %
0260 end
0261 %
0262 % transfer current to hot tail as RF is switched off (TODO : switched down...)
0263 %
0264 cons1t.plh = 0;
0265 cons1t.pecrh = 0;
0266 %
0267 RUNAWAY_EXP.temps = [-1,0];
0268 RUNAWAY_EXP.x = xli;
0269 RUNAWAY_EXP.jrun = RUNAWAY_EXP.jrun + LH_SHAPE_EXP.jlh + ECCD_SHAPE_EXP.jeccd;
0270 RUNAWAY_EXP.irun = [1;1]*(zs.irun + zs.ilh + zs.ieccd);
0271 RUNAWAY_EXP.iohm = [1;1]*NaN;% for first order nonlinear corrections - not used yet
0272 RUNAWAY_EXP.ip = [1;1]*cons1t.ip;
0273 %
0274 disp('Loading external runaway data')
0275 setappdata(0,'RUNAWAY_EXP',RUNAWAY_EXP)
0276 %
0277 [zs,profil,z0dstruct] = zerodevolution(z0dstruct,option,0,cons1t,geo1t,[],sepa1t,[]);
0278 %
0279 % loop on time evolution
0280 %
0281 % time grid (normalized to collision time) -> you can specify :
0282 % - linear arbitrary grid with tn array                 : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0
0283 % - linear arbitrary grid with dtn array                : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10);
0284 % - linear grid with tnmax and grid step                : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000;
0285 % - linear grid with tnmax and number of grid steps     : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i;
0286 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10;
0287 %
0288 if opt.tnmin == 0,%linear time grid
0289     nit = opt.tn/opt.dtn;
0290     tn = linspace(opt.dtn,opt.tn,nit);
0291 end
0292 if opt.tnmin > 0,%log time grid
0293     nit = opt.dtn;
0294     tn = logspace(log10(opt.tnmin),log10(opt.tn),opt.dtn);
0295 end
0296 %
0297 for it = 1:nit,  
0298     %
0299     cons1t.temps = tn(it)/initstate.mksa.nhu_ref;
0300     %
0301     RUNAWAY_EXP.temps(1,end+1) = cons1t.temps;
0302     RUNAWAY_EXP.ip(end+1,1) = cons1t.ip;
0303     RUNAWAY_EXP.iohm(end+1,1) = NaN;% for first order nonlinear corrections - not used yet
0304     %
0305     % loop on E-field convergence
0306     %
0307     RUNAWAY_EXP.jrun(end+1,:) = RUNAWAY_EXP.jrun(end,:);% previous time is first guess
0308     RUNAWAY_EXP.irun(end+1,1) = RUNAWAY_EXP.irun(end,1);% previous time is first guess
0309     %
0310     for iconv = 1:nconv;
0311         %
0312         disp('Loading external runaway data')
0313         setappdata(0,'RUNAWAY_EXP',RUNAWAY_EXP)
0314         %
0315         [zs,profil,z0dstruct] = zerodevolution(z0dstruct,option,cons1t.temps,cons1t,geo1t,[],sepa1t,[]);
0316         %
0317         %
0318     end
0319     
0320     
0321     
0322     
0323     
0324     
0325     
0326     [Znorm(it),Zcurr(it),ZP0(it),dke_out(it)] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,'','',XXfinit{it},[],loadstructs);
0327     XXfinit{it+1} = dke_out(it).XXf0{1};
0328     %
0329     ohm.epsi = dke_out(it).xepsi{end}*mksa0.Edreicer_ref;% guess is previous time
0330     ohm.xnrem_init = dke_out(it).xnorm_rem{end};
0331     ohm.xnrep_init = dke_out(it).xnorm_rep{end};
0332     %
0333     disp(['Time step ',num2str(it),'/',num2str(nit),' calculated.'])
0334 end
0335 
0336 
0337 
0338 
0339 
0340 
0341 
0342 
0343 
0344 
0345 
0346 [dkepath,equil,dkeparam,dkedisplay,waves] = fprep_runaway(wpr,Zi,opt,varargin{:});
0347 %
0348 [~,~,~,~,~,~,~,mc2] = pc_dke_yp;%Universal physics constants
0349 %
0350 %************************************************************************************************************************************
0351 %
0352 % initial (time-asymptotic) distribution
0353 %
0354 equil.pTe = betath^2*mc2*ones(size(equil.pTe));%initial temperature
0355 equil.pzni = equil.pzni*ne0/equil.pne(1);%initial temperature
0356 equil.pne = ne0*ones(size(equil.pne));%initial temperature
0357 %
0358 dkeparam0 = dkeparam;% parameters for initial distribution calculation
0359 %
0360 dkeparam0.tn = NaN;%time-asymptotic
0361 dkeparam0.dtn = NaN;%time-asymptotic
0362 %
0363 % determine the power for LH drive to reach targer current
0364 %
0365 for it_wave = 1:opt.maxiconv,
0366     %
0367     [~,Zcurr0,~,~,~,~,~,~,~,~,~,mksa0,~] = main_dke_yp(id_simul,dkepath,equil,dkeparam0,dkedisplay,'',waves,'','',[],[]);
0368     %
0369     Jtot_rf = Zcurr0.x_0_fsav*mksa0.j_ref;
0370     %
0371     if abs(Jtot_rf - jlh)/jlh <= 1e-3;%opt.tol,
0372         disp(['convergence obtained for RF current after ',num2str(it_wave),' iterations']);
0373         break
0374     elseif it_wave == opt.maxiconv,
0375         disp(['WARNING : convergence not obtained for RF current after ',num2str(it_wave),' iterations']);
0376     else
0377         waves{1}.rays{1}.P0_2piRp = waves{1}.rays{1}.P0_2piRp*jlh/Jtot_rf;
0378     end
0379     %
0380 end
0381 %
0382 dkeparam.Jtot_ohm = j0;% current density
0383 %
0384 dkeparam.nit_ohm = opt.maxiconv;
0385 dkeparam.tol_ohm = opt.tol;
0386 dkeparam.efield_ohm = opt.efield0;
0387 %
0388 dkeparam.opt_struct = 0;%to keep all fields in output structures for loadstruct
0389 %
0390 % determine the initial time asymptotic distribution function to have j0 with ohmic heating
0391 %
0392 dkeparam0 = dkeparam;% parameters for initial distribution calculation
0393 %
0394 dkeparam0.tn = NaN;%time-asymptotic
0395 dkeparam0.dtn = NaN;%time-asymptotic
0396 dkeparam0.opt_struct = 2;%to keep all fields in output structures for loadstruct
0397 %
0398 [Znorm0,Zcurr0,ZP00,dke_out0,radialDKE0,equilDKE0,momentumDKE0,gridDKE0,Zmomcoef0,Zbouncecoef0,Zmripple0,mksa0,XXsinksource0] = ...
0399     main_dke_yp(id_simul,dkepath,equil,dkeparam0,dkedisplay,'','','','',[],[]);
0400 %
0401 dmhu = momentumDKE0.dmhu;
0402 mhu = momentumDKE0.mhu;
0403 dpn = momentumDKE0.dpn;
0404 pn = momentumDKE0.pn;
0405 pn2 = gridDKE0.pn2;
0406 gamma = Zmomcoef0.gamma;
0407 %
0408 if opt.loadstructs,
0409     loadstructs.radialgrid_dke_jd.radialDKE = radialDKE0;
0410     %
0411 %     loadstructs.equilibrium_jd.equilDKE = equilDKE0;
0412     %
0413     loadstructs.momentumgrid_dke_jd.momentumDKE = momentumDKE0;
0414     loadstructs.matrix3Dgridbuilder_dke_yp.gridDKE = gridDKE0;
0415     loadstructs.momentumcoefbuilder_dke_yp.Zmomcoef = Zmomcoef0;
0416     loadstructs.bouncecoefbuilder_dke_yp.Zbouncecoef = Zbouncecoef0;
0417     loadstructs.mripplecoefbuilder_dke_yp.Zmripple = Zmripple0;
0418     loadstructs.mksacoefbuilder_dke_yp.mksa = mksa0;
0419     loadstructs.sinksourcecoeffbuilder_dke_yp.XXsinksource = XXsinksource0;
0420     %
0421     loadstructs.rad_dke_yp.ZXXD_a = dke_out0.ZXXD_a;
0422     loadstructs.rad_dke_yp.ZXXF_a = dke_out0.ZXXF_a;
0423     %
0424 %     loadstructs.coll_dke_jd.ZXXD_c = dke_out0.ZXXD_c;
0425 %     loadstructs.coll_dke_jd.ZXXF_c = dke_out0.ZXXF_c;
0426 %     loadstructs.coll_dke_jd.ZXXD_c_tp = NaN;%dke_out.ZXXD_c_tp;
0427 %     loadstructs.coll_dke_jd.ZXXF_c_tp = NaN;%dke_out.ZXXF_c_tp;
0428 %     loadstructs.coll_dke_jd.XXfM = dke_out0.XXfM;
0429 %     loadstructs.coll_dke_jd.XXILor = dke_out0.XXILor;
0430     %
0431     loadstructs.efield_dke_jd.ZXXD_e = dke_out0.ZXXD_e;
0432     loadstructs.efield_dke_jd.ZXXF_e = dke_out0.ZXXF_e;
0433     loadstructs.efield_dke_jd.ZXXD_e_tp = NaN;%dke_out.ZXXD_e_tp;
0434     loadstructs.efield_dke_jd.ZXXF_e_tp = NaN;%dke_out.ZXXF_e_tp;
0435     loadstructs.efield_dke_jd.xepsi_init = dke_out0.xepsi_init;
0436     loadstructs.efield_dke_jd.xEfield_validity = dke_out0.xEfield_validity;
0437     %
0438 %     loadstructs.rfwave_dke_jd.waveparam = dke_out0.waveparam;
0439 %     loadstructs.rfwave_dke_jd.xyprop_dke = dke_out0.xyprop_dke;
0440 %     loadstructs.rfwave_dke_jd.xyP0_2piRp_mod = dke_out0.xyP0_2piRp_mod;
0441 %     loadstructs.rfwave_dke_jd.yb = dke_out0.yb;
0442 %     loadstructs.rfwave_dke_jd.yP0_2piRp = dke_out0.yP0_2piRp;
0443 %     loadstructs.rfwave_dke_jd.xys = dke_out0.xys;
0444     %
0445 %     loadstructs.loop_rfdiff_dke_jd.ZXYD_rf = dke_out0.ZXYD_rf;
0446 %     loadstructs.loop_rfdiff_dke_jd.ZXYF_rf = dke_out0.ZXYF_rf;
0447 %     loadstructs.loop_rfdiff_dke_jd.ZXYD_rf_tp = NaN;%dke_out.ZXYD_rf_tp;
0448 %     loadstructs.loop_rfdiff_dke_jd.ZXYF_rf_tp = NaN;%dke_out.ZXYF_rf_tp;
0449 %     loadstructs.loop_rfdiff_dke_jd.gridindex_rf = dke_out0.gridindex_rf;
0450     %
0451     % loadstructs.finit_dke_yp.XXfinit = dke_out0.XXfinit;
0452     %
0453 else
0454     %
0455     loadstructs = [];
0456     %
0457 end
0458 %
0459 I_tot = Zcurr0.I_tot;
0460 %
0461 XXfinit = dke_out0.XXf0;
0462 %
0463 ohm = struct;
0464 ohm.id = 'evolution';
0465 ohm.rho = gridDKE0.xrho;
0466 ohm.epsi = dke_out0.xepsi{end}*mksa0.Edreicer_ref;% guess is previous time
0467 ohm.xnrem_init = dke_out0.xnorm_rem{end};
0468 ohm.xnrep_init = dke_out0.xnorm_rep{end};
0469 %
0470 %************************************************************************************************************************************
0471 %
0472 % temperature evolution and time dependent distribution
0473 %
0474 % time grid
0475 %
0476 nit = opt.nit;
0477 %
0478 % time grid (normalized to collision time) -> you can specify :
0479 % - linear arbitrary grid with tn array                 : opt.tnmin = 0 ;opt.tn = 1000:1000:10000,dtn = NaN;opt.tnmin = 0
0480 % - linear arbitrary grid with dtn array                : opt.tnmin = 0 ;opt.tn = NaN;dtn = 1000*ones(1,10);
0481 % - linear grid with tnmax and grid step                : opt.tnmin = 0 ;opt.tn = 10000;dtn = 1000;
0482 % - linear grid with tnmax and number of grid steps     : opt.tnmin = 0 ;opt.tn = 10000;dtn = 10i;
0483 % - log grid with tnmin, tnmax and number of grid steps : opt.tnmin = 1 ;opt.tn = 10000;dtn = 10;
0484 %
0485 dkeparam.tn = NaN;
0486 dkeparam.dtn = opt.tnmax/nit;
0487 tn = linspace(dkeparam.dtn,opt.tnmax,nit);
0488 %
0489 for it = 1:nit,    
0490     [Znorm(it),Zcurr(it),ZP0(it),dke_out(it)] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,'','',XXfinit{it},[],loadstructs);
0491     XXfinit{it+1} = dke_out(it).XXf0{1};
0492     %
0493     ohm.epsi = dke_out(it).xepsi{end}*mksa0.Edreicer_ref;% guess is previous time
0494     ohm.xnrem_init = dke_out(it).xnorm_rem{end};
0495     ohm.xnrep_init = dke_out(it).xnorm_rep{end};
0496     %
0497     disp(['Time step ',num2str(it),'/',num2str(nit),' calculated.'])
0498 end
0499 %
0500 % calculation of the Spitzer current
0501 %
0502 dkeparam.avalanche_mode = 0;
0503 dkeparam.dtn = NaN;%time-asymptotic
0504 dkeparam.nit_ohm = 1;%no field adjustment
0505 %dkeparam = rmfield(dkeparam,{'Te_ref','ne_ref'});%TODO : pb when we keep the old grid...
0506 %
0507 ohm = struct;
0508 ohm.id = 'linear';
0509 ohm.rho = gridDKE0.xrho;
0510 ohm.epsi = 1e-6;% linear limit
0511 %
0512 if opt.loadstructs,
0513     loadstructs = rmfield(loadstructs,'efield_dke_jd');
0514 end
0515 %
0516 [~,Zcurr_lin0] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,'','','',[],[],loadstructs);
0517 %
0518 pfM = cell(1,nit+1);
0519 pf0 = cell(1,nit+1);
0520 %
0521 pfM{1} = (1/2)*integral_dke_jd(dmhu,dke_out0.XXfM,2);%mhu-averaged distribution
0522 pf0{1} = (1/2)*integral_dke_jd(dmhu,dke_out0.XXf0{1},2);%mhu-averaged distribution
0523 %
0524 for it = 1:nit,
0525     pfM{it+1} = (1/2)*integral_dke_jd(dmhu,dke_out(it).XXfM,2);%mhu-averaged distribution
0526     pf0{it+1} = (1/2)*integral_dke_jd(dmhu,dke_out(it).XXf0{1},2);%mhu-averaged distribution
0527 end
0528 %
0529 J.tot = NaN(1,nit+1);% total current density
0530 J.lin = NaN(1,nit+1);% Spitzer current density
0531 J.rei = NaN(1,nit+1);% Internal RE current density
0532 J.reep = NaN(1,nit+1);% External RE current density
0533 J.reem = NaN(1,nit+1);% External RE current density
0534 %
0535 Efield = NaN(1,nit+1);%E/E_c
0536 %
0537 Efield(1) = dke_out0.xepsi{end}/mksa0.betath_ref^2;
0538 %
0539 J.tot(1) = dke_out0.xcurr{end};% total current density
0540 J.lin(1) = Zcurr_lin0.x_0_fsav*(dke_out0.xepsi{end}.*mksa0.Edreicer_ref)/ohm.epsi;% Spitzer current density
0541 J.rei(1) = Zcurr0.x_Rfrac_0_fsav;% Internal RE current density
0542 J.reep(1) = dke_out0.xcurr_rep{end};% External RE current density
0543 J.reem(1) = dke_out0.xcurr_rem{end};% External RE current density
0544 %
0545 for it = 1:nit,
0546     Efield(it+1) = dke_out(it).xepsi{1}/mksa0.betath_ref^2;
0547     %
0548     J.tot(it+1) = dke_out(it).xcurr{end};% total current density
0549     %
0550     % maybe we must multiply Zcurr_lin(it) by the bulk density xnorm0 or ((1-delta)*xnorm0 ?
0551     J.lin(it+1) = Zcurr_lin0.x_0_fsav*(dke_out(it).xepsi{end}*mksa0.Edreicer_ref)/ohm.epsi;% Spitzer current density
0552     J.rei(it+1) = Zcurr(it).x_Rfrac_0_fsav;% Internal RE current density
0553     J.reep(it+1) = dke_out(it).xcurr_rep{end};% External RE current density
0554     J.reem(it+1) = dke_out(it).xcurr_rem{end};% External RE current density
0555 end
0556 %
0557 J.hot = J.tot - J.lin - J.rei - J.reep - J.reem;% Hot tail current density
0558 %
0559 J.j0 = j0/mksa0.j_ref;   
0560 %
0561 [np,nmhu] = size(dke_out0.XXf0{1});
0562 %
0563 savedata = struct;
0564 %
0565 savedata.t = [0,tn];
0566 savedata.betath = betath;
0567 savedata.pn = pn;
0568 savedata.dpn = dpn;
0569 savedata.xi = mhu;
0570 savedata.dxi = dmhu;
0571 %
0572 savedata.f = NaN(np,nmhu,nit+1);
0573 savedata.fM = NaN(np,nmhu,nit+1);
0574 savedata.f(:,:,1) = dke_out0.XXf0{1}/Znorm0.x_0_fsav;
0575 savedata.fM = dke_out0.XXfM;
0576 for it = 1:nit,
0577     savedata.f(:,:,it+1) = dke_out(it).XXf0{1}/Znorm(it).x_0_fsav;
0578     savedata.fM(:,:,it+1) = dke_out(it).XXfM;
0579 end
0580 %
0581 savedata.info.t = 'time grid, normalized to electron thermal collision time 4*pi*e0^2*m^2*vT0^3/(ne*e^4)';
0582 savedata.info.betath = 'ratio of the reference thermal velocity VT0=sqrt(T0/m) to the speed of light c. Also equal to sqrt(T0/mc2)';
0583 savedata.info.pn = 'momentum grid, normalized to the reference thermal momentum pT0 = m*vT0';
0584 savedata.info.dpn = 'momentum grid steps';
0585 savedata.info.xi = 'pitch angle grid, xi = p||/p';
0586 savedata.info.dxi = 'pitch angle grid steps';
0587 savedata.info.rbetath = 'temperature evolution prescribed through field particles in the collision operator: rbetath = vT/vT0';
0588 savedata.info.f = 'electron distribution function as a function of (pn,xi,t), normalized to the electron density';
0589 savedata.info.fM = 'Maxwellian distribution function as a function of (pn,xi,t), normalized to the electron density';
0590 %
frundke_metis

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
