0001 function [lukestructs] = run_lukert(lukestructs,dkepath,timedep)
0002
0003
0004
0005
0006
0007 if nargin < 2
0008 dkepath = load_structures_yp('dkepath','','');
0009 end
0010
0011 if nargin > 2 && timedep && iscell(lukestructs),
0012
0013 nt = length(lukestructs);
0014
0015 for it = 1:nt,
0016
0017 if ~isfield(lukestructs{it}.dkeparam,'opt_flegendre') || lukestructs{it}.dkeparam.opt_flegendre == 0,
0018 disp('Note : option flegendre must be activated for time-dependent LUKE calculations. Enforcing.')
0019 lukestructs{it}.dkeparam.opt_flegendre = 1;
0020 end
0021
0022 if ~isfield(lukestructs{it}.opt,'fields')
0023 lukestructs{it}.opt.fields = {};
0024 end
0025
0026 selectedfields = {'dke_out.XXf0_interp','mksa.Te_ref','mksa.ne_ref','mksa.nhu_ref','radialDKE.xrho_S_dke'};
0027
0028 for ifield = 1:length(selectedfields)
0029 if ~any(strcmp(lukestructs{it}.opt.fields,selectedfields{ifield}))
0030 lukestructs{it}.opt.fields{end+1} = selectedfields{ifield};
0031 end
0032 end
0033
0034 lukestructs{it} = run_lukert(lukestructs{it},dkepath);
0035
0036 if it == nt,
0037 return
0038 end
0039
0040
0041
0042 lukestructs{it+1}.XXf0 = lukestructs{it}.output.dke_out.XXf0_interp{end};
0043
0044
0045
0046
0047
0048
0049 lukestructs{it+1}.dkeparam.Te_ref = lukestructs{it}.output.mksa.Te_ref;
0050 lukestructs{it+1}.dkeparam.ne_ref = lukestructs{it}.output.mksa.ne_ref;
0051 lukestructs{it+1}.dkeparam.rho_S = lukestructs{it}.output.radialDKE.xrho_S_dke;
0052 lukestructs{it+1}.dkeparam.tn = (lukestructs{it+1}.equil.shotime - lukestructs{it}.equil.shotime)*lukestructs{it}.output.mksa.nhu_ref;
0053 lukestructs{it+1}.dkeparam.dtn = 1i;
0054
0055 end
0056
0057 else
0058
0059 cellmode = 1;
0060
0061 if ~iscell(lukestructs) && (any(isfield(lukestructs,{'job','output','err'})) ...
0062 || (isfield(dkepath,'clustermode') && isfield(dkepath.clustermode,'run_lukert'))),
0063
0064 cellmode = 0;
0065 lukestructs = {lukestructs};
0066
0067 end
0068
0069 if iscell(lukestructs),
0070
0071 lukestructsize = size(lukestructs);
0072
0073 lukestructs = lukestructs(:);
0074
0075 njobs = length(lukestructs);
0076
0077 luke_flag = zeros(1,njobs);
0078 statusstr = {'Ready to start','Running','Completed','Failed','Unknown'};
0079
0080 disp(' ')
0081 disp(' -------------------------------------------- ')
0082 disp(' ')
0083
0084 for ijob = 1:njobs,
0085
0086 if isempty(lukestructs{ijob}) && luke_flag(1) == 1,
0087 luke_flag(ijob) = -1;
0088 elseif isfield(lukestructs{ijob},'job'),
0089 luke_flag(ijob) = 1;
0090 elseif isfield(lukestructs{ijob},'output'),
0091 luke_flag(ijob) = 2;
0092 elseif isfield(lukestructs{ijob},'err'),
0093 luke_flag(ijob) = 3;
0094 elseif ~isfield(lukestructs{ijob},'equil'),
0095 luke_flag(ijob) = 4;
0096 end
0097
0098 end
0099
0100 if njobs > 1 && luke_flag(1) == 1 && all(luke_flag(2:end) == -1),
0101 disp(['Checking the status of ',num2str(njobs),' jobs sent to run_lukert via mdce_remote :'])
0102 disp(' ')
0103 else
0104 disp(['Initial status of the ',num2str(njobs),' job(s) sent to run_lukert :'])
0105 disp(' ')
0106
0107 for ijob = 1:njobs,
0108
0109 if luke_flag(ijob) >= 0,
0110 disp([' - Job # ',num2str(ijob),' : ',statusstr{1 + luke_flag(ijob)}])
0111 end
0112 end
0113 end
0114
0115 disp(' ')
0116 disp(' -------------------------------------------- ')
0117 disp(' ')
0118
0119 if any(luke_flag == 0),
0120
0121 if isfield(dkepath,'remnum'),
0122 dkecluster = clustermode_luke('','run_lukert',dkepath,dkepath.remnum);
0123 else
0124 dkecluster = clustermode_luke('','run_lukert',dkepath);
0125 end
0126
0127 [~,lukestructs(luke_flag == 0)] = mdce_run('run_lukert',{''},1,lukestructs(luke_flag == 0),dkecluster);
0128 end
0129
0130 if any(luke_flag == 1),
0131
0132 [~,lukestructs(abs(luke_flag) == 1)] = mdce_get(lukestructs(abs(luke_flag) == 1));
0133
0134 end
0135
0136 if cellmode == 0;
0137 if ~isreal(dkecluster.scheduler.enforce)
0138 lukestructs{1}.warning = ['Remote calculation failed with remote computer number #',int2str(imag(dkecluster.scheduler.enforce)),', local calculation performed'];
0139 end
0140
0141 lukestructs = lukestructs{1};
0142 else
0143 if ~isreal(dkecluster.scheduler.enforce)
0144 for ii = 1:length(lukestructs)
0145 lukestructs{ii}.warning = ['Remote calculation failed with remote computer number #',int2str(imag(dkecluster.scheduler.enforce)),', local calculation performed'];
0146 end
0147 end
0148
0149 lukestructs = reshape(lukestructs,lukestructsize);
0150 end
0151
0152 return
0153
0154 end
0155 end
0156
0157
0158
0159 equil = lukestructs.equil;
0160 dkeparam_in = lukestructs.dkeparam;
0161 dkedisplay_in = lukestructs.dkedisplay;
0162
0163
0164
0165 if ~isfield(lukestructs,'simul'),
0166 simul_in = '';
0167 else
0168 simul_in = lukestructs.simul;
0169 end
0170
0171 if ~isfield(lukestructs,'wavestructs'),
0172 wavestructs = [];
0173 else
0174 wavestructs = lukestructs.wavestructs;
0175 end
0176
0177 if ~isfield(lukestructs,'waves'),
0178 waves = [];
0179 else
0180 waves = lukestructs.waves;
0181 end
0182
0183 if ~isfield(lukestructs,'ohm'),
0184 ohm = [];
0185 else
0186 ohm = lukestructs.ohm;
0187 end
0188
0189 if ~isfield(lukestructs,'transpfaste'),
0190 transpfaste = [];
0191 else
0192 transpfaste = lukestructs.transpfaste;
0193 end
0194
0195 if ~isfield(lukestructs,'ripple'),
0196 ripple = [];
0197 else
0198 ripple = lukestructs.ripple;
0199 end
0200
0201 if ~isfield(lukestructs,'XXf0'),
0202 XXf0 = [];
0203 else
0204 XXf0 = lukestructs.XXf0;
0205 end
0206
0207 if ~isfield(lukestructs,'XXsinksource'),
0208 XXsinksource_in = [];
0209 else
0210 XXsinksource_in = lukestructs.XXsinksource;
0211 end
0212
0213 if ~isfield(lukestructs,'loadstructs'),
0214 loadstructs = [];
0215 else
0216 loadstructs = lukestructs.loadstructs;
0217 end
0218
0219 if ~isfield(lukestructs,'opt'),
0220 opt.save = 0;
0221 opt.fields = NaN;
0222 opt.backup = 0;
0223 opt.waves = 0;
0224 opt.proc = [];
0225 else
0226 opt = lukestructs.opt;
0227 if ~isfield(opt,'save'),
0228 opt.save = 0;
0229 end
0230 if ~isfield(opt,'fields'),
0231 opt.fields = NaN;
0232 end
0233 if ~isfield(opt,'backup'),
0234 opt.backup = 0;
0235 end
0236 if ~isfield(opt,'waves'),
0237 opt.waves = 0;
0238 end
0239 if ~isfield(opt,'proc'),
0240 opt.proc = [];
0241 end
0242 end
0243
0244 [dummy,matlabversion,matlabrelease] = LUKEversion_jd;
0245
0246 if str2double(matlabversion) < 710,
0247
0248
0249
0250 eval('rand(''seed'',sum(100*clock));');
0251 elseif str2double(matlabversion) < 712,
0252
0253
0254
0255 eval('RandStream.setDefaultStream(RandStream(''mt19937ar'',''seed'',sum(100*clock)));');
0256 else
0257 RandStream.setGlobalStream(RandStream('mt19937ar','seed',sum(100*clock)));
0258 end
0259
0260 if ~isfield(equil,'fluct') || isempty(equil.fluct),
0261 equil.fluct = struct;
0262 end
0263
0264 if ~isstruct(simul_in),
0265 simul.id = simul_in;
0266 simul.locid = '';
0267 simul.path = './';
0268 else
0269 simul = simul_in;
0270 if ~isfield(simul,'id'),
0271 simul.id = '';
0272 end
0273 if ~isfield(simul,'locid'),
0274 simul.locid = '';
0275 end
0276 if ~isfield(simul,'path'),
0277 simul.path = './';
0278 end
0279 end
0280
0281 if isempty(simul.id),
0282 simul.id = make_luke_simulid_jd(equil,ohm,transpfaste,ripple,waves,wavestructs,simul.locid);
0283 end
0284
0285 if ~exist(simul.path,'dir'),
0286 if ~isempty(simul.path),
0287 mkdir(simul.path);
0288 end
0289 end
0290
0291 if opt.waves == 1 && ~exist([simul.path,'wave_results'],'dir'),
0292 mkdir([simul.path,'wave_results']);
0293 end
0294
0295 if opt.backup < 0 && ~exist([simul.path,'backup_',simul.id,'.mat'],'file'),
0296 opt.backup = 1;
0297 end
0298
0299 disp(' ');
0300 disp(['Simulation id : ',simul.id]);
0301 disp(['Simulation path : ',simul.path]);
0302 disp(' ');
0303
0304 itnmin = 1;
0305 wavetype = 0;
0306
0307 if ~isempty(waves) && isstruct(waves{1}),
0308 for iw = 1:length(waves),
0309 if strcmp(waves{iw}.equil_id,equil.id) == 0,
0310 error('The wave structure was build with an inconsistent equilibrium')
0311 end
0312 end
0313 end
0314
0315 f0struct = [];
0316
0317 if isempty(dkeparam_in.psin_S) && isempty(dkeparam_in.rho_S),
0318 dkeparam_in.rt_mode = 0;
0319 disp('WARNING: dkeparam.rt_mode = 0 has been enforced, since both psin_S and rho_S are empty.');
0320 end
0321
0322 if ~isfield(dkeparam_in,'tn') || (isempty(dkeparam_in.psin_S) && isempty(dkeparam_in.rho_S)),
0323 dkeparam_in.tn = NaN;
0324 disp('WARNING: dkeparam.tn = NaN has been set, since tn field is not existing in dkeparam structure.');
0325 end
0326
0327 if ~isfield(dkeparam_in,'Nfluct'),
0328 dkeparam_in.Nfluct = 1;
0329 end
0330
0331
0332
0333 flag_LH = 0;
0334
0335 for iw = 1:length(wavestructs),
0336 if strcmp(wavestructs{iw}.launch.type,'LH'),
0337 flag_LH = 1;
0338 end
0339 end
0340
0341 if flag_LH == 0 && isfield(equil.fluct,'npar0_lh'),
0342 disp('Perturbation of the launched npar for the LH wave only; dismissed.')
0343 equil.fluct = rmfield(equil.fluct,'npar0_lh');
0344 end
0345
0346 mksa = buildlukestruct_yp('mksa',equil,dkeparam_in,dkedisplay_in);
0347
0348 if isfield(dkeparam_in,'tn_mode') && dkeparam_in.tn_mode == 1,
0349
0350 if ~isscalar(dkeparam_in.tn) || (~isnan(dkeparam_in.tn) && imag(dkeparam_in.tn) == 0),
0351 dkeparam_in.tn = mksa.nhu_ref*dkeparam_in.tn;
0352 end
0353 if ~isscalar(dkeparam_in.dtn) || (~isnan(dkeparam_in.dtn) && imag(dkeparam_in.dtn) == 0),
0354 dkeparam_in.dtn = mksa.nhu_ref*dkeparam_in.dtn;
0355 end
0356
0357 end
0358
0359
0360
0361 [dtnmin_fluct,dtns_fluct] = calc_dtn_fluct_jd(equil.fluct,wavestructs,mksa);
0362
0363 dtn_phase = dkeparam_in.Nfluct*dtnmin_fluct;
0364
0365 if isfield(equil.fluct,'itnmax'),
0366 itnmax = equil.fluct.itnmax;
0367 tn_in = itnmax*dtn_phase;
0368 dkeparam_in.dtn = 1i;
0369 elseif isscalar(dkeparam_in.tn) && imag(dkeparam_in.tn) > 0,
0370 itnmax = imag(dkeparam_in.tn);
0371 tn_in = itnmax*dtn_phase;
0372 elseif isscalar(dkeparam_in.tn) && ~isnan(dkeparam_in.tn),
0373 itnmax = max(1,ceil(dkeparam_in.tn/dtn_phase));
0374 tn_in = dkeparam_in.tn;
0375 else
0376 itnmax = 1;
0377 tn_in = dkeparam_in.tn;
0378 end
0379
0380 dkeparam_in = rmfield(dkeparam_in,'tn');
0381
0382 byPnonabs = NaN;
0383
0384 if itnmax > 1 && ~isempty(waves),
0385 twaves = waves;
0386 clear waves
0387 else
0388 twaves = [];
0389 end
0390
0391 if opt.backup >= 0;
0392
0393 dkeparam = dkeparam_in;
0394 dkedisplay = dkedisplay_in;
0395 tn = tn_in;
0396
0397 calctime = struct;
0398
0399 if ~isinf(dtnmin_fluct),
0400 xJ = [];
0401 xP = [];
0402 Itot = [];
0403 Ptot = [];
0404 end
0405
0406
0407
0408 if isstruct(equil.fluct) && ~isempty(fieldnames(equil.fluct)),
0409 equil.fluct = fluctphase_yp(equil.fluct);
0410 end
0411
0412 if ~isempty(wavestructs),
0413 for iws = 1:length(wavestructs),
0414 [wavestructs{iws}.equil] = equil;
0415 if isfield(wavestructs{iws},'fitparam'),
0416 [wavestructs{iws}.equil_fit] = fitequil_yp(equil,wavestructs{iws}.fitparam.mode_equil,wavestructs{iws}.fitparam.method,wavestructs{iws}.fitparam.ngridresample,wavestructs{iws}.fitparam.nharm);
0417 else
0418 [wavestructs{iws}.equil_fit] = fitequil_yp(equil,wavestructs{iws}.fitequilparam.mode_equil,wavestructs{iws}.fitequilparam.method,wavestructs{iws}.fitequilparam.ngridresample,wavestructs{iws}.fitequilparam.nharm);
0419 end
0420
0421 if isstruct(equil.fluct) && ~isempty(fieldnames(equil.fluct)),
0422 if isfield(wavestructs{iws},'fitparam'),
0423 [wavestructs{iws}.equil_fit.fluct_fit] = fitfluct_yp(equil.fluct,wavestructs{iws}.fitparam.mode_equil,equil.fluct.fitparam.method,equil.fluct.fitparam.ngridresample,equil.fluct.fitparam.nharm);
0424 else
0425 if isfield(equil.fluct,'fitparam'),
0426 [wavestructs{iws}.equil_fit.fluct_fit] = fitfluct_yp(equil.fluct,wavestructs{iws}.fitequilparam.mode_equil,equil.fluct.fitparam.method,equil.fluct.fitparam.ngridresample,equil.fluct.fitparam.nharm);
0427 else
0428 [wavestructs{iws}.equil_fit.fluct_fit] = fitfluct_yp(equil.fluct);
0429 end
0430 end
0431 end
0432 end
0433 end
0434 else
0435
0436 disp('------------------ RE-LOAD MODE ------------------');
0437
0438 load([simul.path,'backup_',simul.id,'.mat']);
0439
0440 if opt.backup == -1,
0441
0442 dkeparam = imod_struct_jd(dkeparam,'dkeparam',1,dkedisplay.display_mode);
0443 dkedisplay = imod_struct_jd(dkedisplay,'dkedisplay',1,dkedisplay.display_mode);
0444 opt = imod_struct_jd(opt,'opt',1,dkedisplay.display_mode);
0445
0446 for ifluct = 1:size(wavestructs,1),
0447 for iw = 1:size(wavestructs,2),
0448
0449 wavestructs{ifluct,iw} = imod_struct_jd(wavestructs{ifluct,iw},['wavestructs{',num2str(ifluct),',',num2str(iw),'}'],1,dkedisplay.display_mode);
0450 end
0451 end
0452
0453 elseif opt.backup == -2,
0454 dkedisplay.display_mode = dkedisplay_in.display_mode;
0455 if ~isfield(dkedisplay_in,'display_time_mode');
0456 dkedisplay.display_time_mode = 0;
0457 else
0458 dkedisplay.display_time_mode = dkedisplay_in.display_time_mode;
0459 end
0460 end
0461
0462 if itn == itnmax,
0463 disp(['The maximum number of iterations ',int2str(itnmax),' has already been reached.']);
0464 else
0465 disp(['The number of iterations already done is ',int2str(itn),'/',int2str(itnmax),'.']);
0466 end
0467
0468 if tn < tn_in,
0469 if opt.backup == -1,
0470 overwrite_tn = input_dke_yp(['Do you want to overwrite the total simulation time from ',int2str(tn),' to ',int2str(tn_in),' (y/n) ?'],'n',{'y','n'});
0471 else
0472 overwrite_tn = 'y';
0473 end
0474
0475 if overwrite_tn == 'y',
0476
0477 disp(['- Old tn = ',int2str(tn)]);
0478 disp(['- New tn = ',int2str(tn_in)]);
0479 tn = tn_in;
0480 itnmax = max([1,ceil(tn/dtn_phase)]);
0481 end
0482
0483 end
0484
0485 opt.backup = 1;
0486
0487 disp('---------------------------------')
0488 disp('Iteration status')
0489 disp('---------------------------------')
0490 disp([' - tn = ',int2str(tn)])
0491 disp([' - number of iterations to be done = ',int2str(itnmax - itn)])
0492
0493 if itn < itnmax,
0494 itnmin = itn + 1;
0495 else
0496 itnmin = itnmax + 1;
0497 end
0498 end
0499
0500 for itn = itnmin:itnmax,
0501
0502 if itn < itnmax,
0503 dkeparam.tn = dtn_phase;
0504 elseif isscalar(tn) && ~isinf(dtn_phase),
0505 dkeparam.tn = tn - dtn_phase*(itnmax - 1);
0506 else
0507 dkeparam.tn = tn;
0508 end
0509
0510
0511
0512 time_waves = tic;
0513
0514 if ~isempty(wavestructs),
0515 [waves,wavestructs,byisel] = wave_distribution_yp(dkepath,dkeparam,wavestructs,f0struct,dtnmin_fluct/mksa.nhu_ref,byPnonabs);
0516 elseif ~isempty(twaves),
0517 [waves,byisel] = wave_power_jd(dkeparam,twaves{itn},byPnonabs);
0518 end
0519
0520 calctime(itn).waves = toc(time_waves);
0521
0522
0523
0524
0525 if itn > 1,
0526 XXf0 = XXf0_t{itn-1};
0527 end
0528
0529 time_luke = tic;
0530
0531 if itn == 1,
0532 [Znorm_loc,Zcurr_loc,ZP0_loc,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = loop_run_lukert(simul,waves,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,XXf0,XXsinksource_in,loadstructs);
0533
0534 equilDKE.fluct = equil.fluct;
0535
0536 loadstructs.radialgrid_dke_jd.radialDKE = radialDKE;
0537 loadstructs.equilibrium_jd.equilDKE = equilDKE;
0538 loadstructs.momentumgrid_dke_jd.momentumDKE = momentumDKE;
0539 loadstructs.matrix3Dgridbuilder_dke_yp.gridDKE = gridDKE;
0540 loadstructs.momentumcoefbuilder_dke_yp.Zmomcoef = Zmomcoef;
0541 loadstructs.bouncecoefbuilder_dke_yp.Zbouncecoef = Zbouncecoef;
0542 loadstructs.mripplecoefbuilder_dke_yp.Zmripple = Zmripple;
0543 loadstructs.mksacoefbuilder_dke_yp.mksa = mksa;
0544 loadstructs.sinksourcecoeffbuilder_dke_yp.XXsinksource = XXsinksource;
0545
0546 loadstructs.rad_dke_yp.ZXXD_a = dke_out.ZXXD_a;
0547 loadstructs.rad_dke_yp.ZXXF_a = dke_out.ZXXF_a;
0548
0549 loadstructs.coll_dke_jd.ZXXD_c = dke_out.ZXXD_c;
0550 loadstructs.coll_dke_jd.ZXXF_c = dke_out.ZXXF_c;
0551 loadstructs.coll_dke_jd.ZXXD_c_tp = NaN;
0552 loadstructs.coll_dke_jd.ZXXF_c_tp = NaN;
0553 loadstructs.coll_dke_jd.XXfM = dke_out.XXfM;
0554 loadstructs.coll_dke_jd.XXILor = dke_out.XXILor;
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578 else
0579 [Znorm_loc,Zcurr_loc,ZP0_loc,dke_out] = loop_run_lukert(simul,waves,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,XXf0,XXsinksource_in,loadstructs);
0580 end
0581
0582 calctime(itn).luke = toc(time_luke);
0583
0584 Ptot_0 = 0;
0585 Ptot_inj = 0;
0586 Pnonabs = 0;
0587
0588 if isfield(dkeparam,'enforcepower') && dkeparam.enforcepower,
0589 byPnonabs = cell(size(wavestructs));
0590
0591 if ~isinf(dtn_phase),
0592
0593 nit_rf = length(dke_out.residue_rf{end});
0594
0595 for ifluct = 1:size(wavestructs,1),
0596 for iw = 1:size(wavestructs,2),
0597
0598
0599
0600 if isfield(wavestructs{ifluct,iw}.launch,'bPlhtot'),
0601 yPtot_0 = wavestructs{ifluct,iw}.launch.bPlhtot_in;
0602 yPtot_inj = wavestructs{ifluct,iw}.launch.bPlhtot;
0603 elseif isfield(wavestructs{ifluct,iw}.launch,'yP_L'),
0604 yPtot_0 = wavestructs{ifluct,iw}.launch.yP_L_in;
0605 yPtot_inj = wavestructs{ifluct,iw}.launch.yP_L;
0606 end
0607
0608 Ptot_0 = Ptot_0 + sum(yPtot_0)/1e6;
0609 Ptot_inj = Ptot_inj + sum(yPtot_inj)/1e6;
0610 ny = length(yPtot_0);
0611
0612
0613
0614 byPnonabs{ifluct,iw} = yPtot_inj;
0615 for iy = 1:ny,
0616 for iray = 1:length(byisel{ifluct,iw}{iy}),
0617 Pabs_loc = 2*pi*equil.Rp*(dke_out.yP0_2piRp(byisel{ifluct,iw}{iy}(iray)) - dke_out.zP_0_2piRp_mod{nit_rf,end}{byisel{ifluct,iw}{iy}(iray)}(end));
0618 byPnonabs{ifluct,iw}(iy) = byPnonabs{ifluct,iw}(iy) - Pabs_loc;
0619 end
0620 end
0621
0622 Pnonabs = Pnonabs + sum(byPnonabs{ifluct,iw})/1e6;
0623
0624 end
0625 end
0626
0627 disp(['---> Injected rf power @ itn = 1 : ',num2str(Ptot_0),' (MW)']);
0628 disp(['---> Injected rf power @ itn = ',int2str(itn),' : ',num2str(Ptot_inj),' (MW)']);
0629 disp(['---> Power not absorbed @ itn = ',int2str(itn),' : ',num2str(Pnonabs),' (MW)']);
0630 end
0631 else
0632 for iw = 1:length(waves),
0633 for ib = 1:length(waves{iw}.rays),
0634 Ptot_inj = Ptot_inj + 2*pi*equil.Rp*sum(waves{iw}.rays{ib}.P0_2piRp)/1e6;
0635 end
0636 end
0637 end
0638
0639 if itnmax > 1,
0640
0641 xrho = loadstructs.equilibrium_jd.equilDKE.xrho(:);
0642 xdA = loadstructs.equilibrium_jd.equilDKE.xdA_dke(:).';
0643 if isfield(loadstructs.equilibrium_jd.equilDKE,'xdV_dke'),
0644 xdV = loadstructs.equilibrium_jd.equilDKE.xdV_dke(:).';
0645 else
0646 xdV = loadstructs.equilibrium_jd.equilDKE.xdV_2piRp_dke(:).'*2*pi*loadstructs.equilibrium_jd.equilDKE.Rp;
0647 end
0648
0649 xtn(itn) = min([dtn_phase*itn,tn]);
0650
0651 xJ(itn,:) = Zcurr_loc(end).x_0_fsav*loadstructs.mksacoefbuilder_dke_yp.mksa.j_ref;
0652 xP(itn,:) = ZP0_loc(end).x_rf_fsav.'*loadstructs.mksacoefbuilder_dke_yp.mksa.P_ref;
0653
0654 Itot(itn) = sum(xJ(itn,:).*xdA,2);
0655 Ptot(itn) = sum(xP(itn,:).*xdV,2);
0656
0657 if ~isfield(dkedisplay,'display_time_mode'), dkedisplay.display_time_mode = 0;end
0658
0659 if dkedisplay.display_time_mode >= 1,
0660
0661 figure(figbyname_yp(['Jrf(r,t) / ',simul.id])),
0662
0663 if itn > 2,
0664 graph1D_jd(xrho,xJ([1:end-2],:),0,0,'\rho','J_{tot} (MA/m^2)',['Iteration = ',int2str(itn),'/',int2str(itnmax),', I_{tot} = ',num2str(Itot(itn)),' (MA)'],NaN,'','','-','none','g',1);
0665 end
0666
0667 if itn > 1,
0668 graph1D_jd(xrho,xJ(end-1,:),0,0,'\rho','J_{tot} (MA/m^2)',['Iteration = ',int2str(itn),'/',int2str(itnmax),', I_{tot} = ',num2str(Itot(itn)),' (MA)'],NaN,'','','-','none','b',1);
0669 end
0670
0671 graph1D_jd(xrho,xJ(end,:),0,0,'\rho','J_{tot} (MA/m^2)',['I_{tot} = ',num2str(Itot(itn)),' (MA)'],NaN,'','','-','none','r',2);
0672
0673 figure(figbyname_yp(['Prf(r,t) / ',simul.id])),
0674
0675 if itn > 2,
0676 graph1D_jd(xrho,xP([1:end-2],:),0,0,'\rho','P_{abs} (MW/m^3)',['Iteration = ',int2str(itn),'/',int2str(itnmax),', P_{abs} = ',num2str(Ptot(itn)),' (MW)'],NaN,'','','-','none','g',1);
0677 end
0678
0679 if itn > 1,
0680 graph1D_jd(xrho,xP(end-1,:),0,0,'\rho','P_{abs} (MW/m^3)',['Iteration = ',int2str(itn),'/',int2str(itnmax),', P_{abs} = ',num2str(Ptot(itn)),' (MW)'],NaN,'','','-','none','b',1);
0681 end
0682
0683 graph1D_jd(xrho,xP(end,:),0,0,'\rho','P_{abs} (MW/m^3)',['P_{abs} = ',num2str(Ptot(itn)),' (MW)'],NaN,'','','-','none','r',2);
0684
0685 figure(figbyname_yp(['Irf(t) / ',simul.id])),
0686
0687 graph1D_jd(xtn,Itot,0,0,'t*\nu','I_{tot} (MA)',['Iteration = ',int2str(itn),'/',int2str(itnmax)],NaN,'','','none','.','c',2);
0688 graph1D_jd(xtn,Itot,0,0,'t*\nu','I_{tot} (MA)',['Iteration = ',int2str(itn),'/',int2str(itnmax)],NaN,'','','-','none','c',2);
0689
0690 figure(figbyname_yp(['Prf(t) / ',simul.id])),
0691
0692 graph1D_jd(xtn,Ptot,0,0,'t*\nu','P_{tot} (MW)',['Iteration = ',int2str(itn),'/',int2str(itnmax)],NaN,'',[0,Ptot_inj*2],'none','.','k',2);drawnow
0693 graph1D_jd(xtn,Ptot,0,0,'t*\nu','P_{tot} (MW)',['Iteration = ',int2str(itn),'/',int2str(itnmax)],NaN,'',[0,Ptot_inj*2],'-','none','k',2);drawnow
0694 graph1D_jd(xtn,Ptot_inj*ones(1,itn),0,0,'t*\nu','P_{tot} (MW)',['Iteration = ',int2str(itn),'/',int2str(itnmax)],NaN,'',[0,Ptot_inj*2],'--','none','k',1);drawnow
0695 else
0696 disp(' ');
0697 disp('------------------------------------------------------------------------------------');
0698 disp(['Simulation ',simul.id]);
0699 disp(['Iteration #',int2str(itn),': Pref = ',num2str(Ptot_inj),' (MW) - Ptot = ',num2str(Ptot(itn)),' (MW) - Itot = ',num2str(Itot(itn)),' (MA)']);
0700 disp('------------------------------------------------------------------------------------');
0701 disp(' ');
0702 end
0703
0704
0705
0706 f0struct = build_f0struct_yp(dke_out,loadstructs.momentumgrid_dke_jd.momentumDKE,loadstructs.equilibrium_jd.equilDKE,loadstructs.mksacoefbuilder_dke_yp.mksa);
0707 end
0708
0709 XXf0_t{itn} = dke_out.XXf0{end};
0710
0711
0712
0713 for iws = 1:length(wavestructs),
0714 if isfield(wavestructs{iws}.launch,'yP_L');
0715 wavetype = wavetype + 1;
0716 end
0717 end
0718
0719
0720
0721 if itn > 1,
0722 Znorm = [Znorm,Znorm_loc];
0723 Zcurr = [Zcurr,Zcurr_loc];
0724 ZP0 = [ZP0,ZP0_loc];
0725
0726 Ztime = [Ztime,Ztime(end) + dke_out.FPtime];
0727 Zitn = [Zitn,length(Ztime)];
0728 else
0729 Znorm = Znorm_loc;
0730 Zcurr = Zcurr_loc;
0731 ZP0 = ZP0_loc;
0732
0733 Ztime = dke_out.FPtime;
0734 Zitn = length(Ztime);
0735 end
0736
0737 if itnmax > 1 && opt.backup > 0,
0738 wavetype = wavetype/length(wavestructs);
0739
0740 if isfield(dkedisplay,'display_time_mode'),
0741 if dkedisplay.display_time_mode >= 2,
0742 rayevolution_yp(wavetype,wavestructs{1}.equil_fit,waves,0,loadstructs.radialgrid_dke_jd.radialDKE.xrho_S_dke,itn);
0743 end
0744 end
0745
0746 if opt.waves == 1,
0747 save([simul.path,'wave_results/waves_',simul.id,'_',int2str(itn),'.mat'],'wavetype','waves','wavestructs','itn');
0748 end
0749
0750 wavetype = 0;
0751
0752
0753
0754 save_str = [simul.path,'backup_',simul.id,'.mat'];
0755 if isfield(dke_out,'waves'),
0756 waves = dke_out.waves;
0757 else
0758 waves = '';
0759 end
0760
0761 eval(['save ',save_str,' tn itn itnmax calctime Ztime Zitn xrho xtn xJ xP Itot Ptot byPnonabs wavestructs waves f0struct XXf0_t loadstructs simul equil dkeparam dkedisplay ohm transpfaste ripple XXsinksource_in Znorm Zcurr ZP0']);
0762 disp('--------------------------------------------------------------------------------------------------');
0763 info_dke_yp(2,['Backup of iteration ',int2str(itn),'/',int2str(itnmax),' done in file ''',save_str,'''.']);
0764 disp('--------------------------------------------------------------------------------------------------');
0765 clear waves
0766
0767 elseif itnmax == 1 && opt.waves == 1
0768
0769 save([simul.path,'wave_results/waves_',simul.id,'.mat'],'waves','wavestructs');
0770
0771 end
0772 end
0773
0774 dke_out.dkeparam = dkeparam;
0775 dke_out.calctime = calctime;
0776
0777 output.Znorm = Znorm;
0778 output.Zcurr = Zcurr;
0779 output.ZP0 = ZP0;
0780
0781 test_legendre = 0;
0782 if isfield(dkeparam,'opt_flegendre') && dkeparam.opt_flegendre == 1,
0783 [dke_out.XXf0_interp,testf0] = interp_legendre_yp(dkeparam,gridDKE,dke_out.XXf0,test_legendre);
0784 [dke_out.XXfM_interp,testfM] = interp_legendre_yp(dkeparam,gridDKE,dke_out.XXfM,test_legendre);
0785 [dke_out.XXfinit_interp,testfinit] = interp_legendre_yp(dkeparam,gridDKE,dke_out.XXfinit,test_legendre);
0786
0787
0788
0789 end
0790
0791 output.dke_out = dke_out;
0792 output.radialDKE = loadstructs.radialgrid_dke_jd.radialDKE;
0793 output.equilDKE = loadstructs.equilibrium_jd.equilDKE;
0794 output.momentumDKE = loadstructs.momentumgrid_dke_jd.momentumDKE;
0795 output.gridDKE = loadstructs.matrix3Dgridbuilder_dke_yp.gridDKE;
0796 output.Zmomcoef = loadstructs.momentumcoefbuilder_dke_yp.Zmomcoef;
0797 output.Zbouncecoef = loadstructs.bouncecoefbuilder_dke_yp.Zbouncecoef;
0798 output.Zmripple = loadstructs.mripplecoefbuilder_dke_yp.Zmripple;
0799 output.mksa = loadstructs.mksacoefbuilder_dke_yp.mksa;
0800 output.XXsinksource = loadstructs.sinksourcecoeffbuilder_dke_yp.XXsinksource;
0801
0802 if abs(opt.waves),
0803 output.dke_out = rmfield(output.dke_out,'waves');
0804 end
0805
0806
0807
0808 if ~isfield(dkeparam,'opt_struct') || dkeparam.opt_struct == 0,
0809 savefields_radialDKE = {'xrho_S_dke','xpsin_f','xpsin_f_dke','xdpsin_f_dke','xpsin_S_dke','r_dke'};
0810 savefields_momentumDKE = {'pn','pnp','pn_S','dpn','mhu','mhu_S','dmhu'};
0811 savefields_gridDKE = {'X1mmhu2','xmhubounce2'};
0812 savefields_Zmomcoef = {'gamma','v','Ec_ref','dEc_ref'};
0813 if isfield(output.dke_out,'xepsi')
0814 if mean(output.dke_out.xepsi{1}) > 0,
0815 savefields_Zbouncecoef = {'xq','xqtilde','xqbar','xqhat','XXheaviside','Xxlambda','XXlambda_b_p1m2p2'};
0816 else
0817 savefields_Zbouncecoef = {'xq','xqtilde','xqbar','xqhat','XXheaviside','Xxlambda'};
0818 end
0819 end
0820 savefields_Zmripple = {'XXlossripple','xdeltaripple','xNq','Xnhuloss','xthetabanana','xthetaloss','xpndetrap','xpnth','xEseuil'};
0821
0822 savefields_ZXXD_rf = {'parpar_ij','pp_ij'};
0823 rmfields_dke_out = {'ZXXS','ZXXF_rf','ZXXD_c','ZXXF_c','ZXXD_e','ZXXF_e','ZXXD_a','ZXXF_a','XXILor'};
0824
0825 output.radialDKE = selectfields_jd(output.radialDKE,savefields_radialDKE);
0826 output.momentumDKE = selectfields_jd(output.momentumDKE,savefields_momentumDKE);
0827 output.gridDKE = selectfields_jd(output.gridDKE,savefields_gridDKE);
0828 output.Zmomcoef = selectfields_jd(output.Zmomcoef,savefields_Zmomcoef);
0829 output.Zbouncecoef = selectfields_jd(output.Zbouncecoef,savefields_Zbouncecoef);
0830 output.Zmripple = selectfields_jd(output.Zmripple,savefields_Zmripple);
0831
0832 output.dke_out.ZXXD_rf = selectfields_jd(output.dke_out.ZXXD_rf,savefields_ZXXD_rf);
0833 output.dke_out = rmfield(output.dke_out,rmfields_dke_out);
0834 end
0835
0836 lukestructs.output = output;
0837
0838 if opt.save == 1,
0839 if isempty(simul.path) || ~strcmp(simul.path(1),filesep),
0840 simul.path = [pwd,filesep,simul.path];
0841 end
0842 savepath = simul.path;
0843 savefile = ['LUKE_RESULTS_',simul.id,'.mat'];
0844 save([savepath,savefile],'lukestructs','-v7.3');
0845 end
0846
0847 if opt.save >= 0 && isfield(opt,'proc') && isstruct(opt.proc),
0848 lukestructs.output.data_proc = proc_luke_jd(output,opt.proc);
0849 if iscell(opt.fields),
0850 opt.fields = [opt.fields,'data_proc'];
0851 else
0852 opt.fields = {'data_proc'};
0853 end
0854 end
0855
0856 if iscell(opt.fields),
0857 lukestructs.output = selectfields_jd(lukestructs.output,opt.fields);
0858 end
0859
0860 if opt.save == 1,
0861 lukestructs.output.savepath = savepath;
0862 lukestructs.output.savefile = savefile;
0863 end
0864
0865
0866