metis2luke

 function for running LUKE at a given time in the METIS profil0d time grid

 by J. Decker et J.-F. Artaud, 02/01/2011

0001 function output = metis2luke(post,it,auto_mode)
0002 %
0003 % function for running LUKE at a given time in the METIS profil0d time grid
0004 %
0005 % by J. Decker et J.-F. Artaud, 02/01/2011
0006 %
0007 output = '';
0008 %
0009 profil0d = post.profil0d;
0010 zerod = post.zerod;
0011 exp0d = post.z0dinput.exp0d;
0012 cons = post.z0dinput.cons;
0013 geo = post.z0dinput.geo;
0014 option = post.z0dinput.option;
0015 tokamak = post.z0dinput.machine;
0016 shotnum = num2str(fix(abs(post.z0dinput.shot)));
0017 %
0018 % clear post
0019 %
0020 % test inputs and set calculation time
0021 %
0022 if nargin < 3,
0023     auto_mode = 0;
0024 end
0025 %
0026 if nargin < 2,
0027     trange = [min(profil0d.temps);max(profil0d.temps)];
0028     t = input_dke_yp(['Please enter the discharge time at which LUKE is to be called (s) [',num2str(trange.'),']'],mean(trange),trange);
0029     it = find(abs(t - profil0d.temps) == min(abs(t - profil0d.temps)),1,'first');
0030     if profil0d.temps(it) ~= t,
0031         t = profil0d.temps(it);
0032         disp(['Closest METIS profiles time selected : t = ',num2str(t),' s.']);
0033     end
0034 else
0035     if ~isscalar(it) || it ~= round(it) || it < 1 || it > length(profil0d.temps),
0036         error('Invalid METIS profiles time index selected')
0037     else
0038         t = profil0d.temps(it);
0039         disp(['METIS profiles time selected : t = ',num2str(t),' s.']);
0040     end
0041 end
0042 %
0043 shotime = num2str(t,'%6.4f');
0044 %
0045 % extract METIS data at computation time
0046 %
0047 profil0d = structindex(profil0d,length(profil0d.temps),it,{'xli'});
0048 %
0049 itzerod = find(abs((zerod.temps - t)/t) < 10*eps,1,'first');
0050 if isempty(itzerod),
0051     error('profil0d time not found in exp0d time grid.')
0052 else
0053     zerod = structindex(zerod,length(zerod.temps),itzerod);
0054 end
0055 %
0056 itexp = find(abs((exp0d.temps - t)/t) < 10*eps,1,'first');
0057 if isempty(itexp),
0058     error('profil0d time not found in exp0d time grid.')
0059 else
0060     exp0d = structindex(exp0d,length(exp0d.temps),itexp);
0061 end
0062 %
0063 itcons = find(abs((cons.temps - t)/t) < 10*eps,1,'first');
0064 if isempty(itcons),
0065     error('profil0d time not found in cons time grid.')
0066 else
0067     ntcons = length(cons.temps);
0068     cons = structindex(cons,ntcons,itcons);
0069     geo = structindex(geo,ntcons,itcons);
0070 end
0071 %
0072 % LUKE parameters structure
0073 %
0074 lukeparam = make_lukeparam_metis2luke(cons,option,auto_mode);
0075 dkeparam = lukeparam.dkeparam;
0076 %
0077 % LUKE display structure
0078 %
0079 dkedisplay = make_display_cronos;
0080 dkedisplay.rho_display  = lukeparam.rho_display;%Approximate radial position at which results are displayed if display.display_mode >= 1 [1,1]. Not used if NaN
0081 dkedisplay.display_mode = lukeparam.display_mode;%Display mode: (0): no, (1): partial, (2): full
0082 %
0083 % LUKE equilibrium structure
0084 %
0085 equil_id = [tokamak,'_',shotnum,'_',shotime];
0086 %
0087 [equil,profil0d_out] = make_equil_metis2luke(equil_id,profil0d,zerod,exp0d,cons,geo,option,lukeparam.equilparam);
0088 %
0089 prho = equil.ptx(:,1).'/equil.ptx(end,1);%rho grid in LUKE definition
0090 %
0091 % LUKE wavestructs structure
0092 %
0093 launchs_lh = make_wavelaunchs_metis2luke_lh(option,lukeparam.lhparam);
0094 launchs_ec = make_wavelaunchs_metis2luke_ec(option,lukeparam.ecparam,equil);
0095 %
0096 launchs = [launchs_lh,launchs_ec];
0097 %
0098 nw = length(launchs);
0099 %
0100 wavestructs = cell(1,nw);
0101 %
0102 dkeparam.w_mask = [];% waves to be accounted for in radial grid optimization: EC only or LH if no EC present
0103 %
0104 simul_wid = '';
0105 %
0106 for iw = 1:nw,
0107     wavestructs{iw} = make_wavert_cronos(launchs{iw});%parameter for ray-tracing are adapted to the selected wave
0108     %
0109     wavestructs{iw}.equil_id = equil.id;
0110     %
0111     if strcmp(launchs{iw}.type,'EC')
0112         %
0113         dkeparam.w_mask = [dkeparam.w_mask,iw];
0114         disp(['wave :',int2str(iw),' type :',launchs{iw}.type])
0115         %
0116     else
0117         disp(['wave :',int2str(iw),' type :',launchs{iw}.type,'; n// = ',num2str(wavestructs{iw}.launch.bNpar0,4),...
0118             '; Plh = ',num2str(wavestructs{iw}.launch.bPlhtot/1e6,4),' MW'])
0119     end
0120     %
0121     simul_wid = [simul_wid,'_',wavestructs{iw}.id];
0122     %
0123 end
0124 if isempty(dkeparam.w_mask)
0125     dkeparam.w_mask = 1:length(launchs);
0126 end
0127 %
0128 % LUKE ohm structure
0129 %
0130 ohm.id = 'EPAR';
0131 ohm.rho = prho;
0132 ohm.epsi = - profil0d_out.epar;%in V/m (the minus sign is for the electron charge)
0133 %
0134 % LUKE radial transport structure
0135 %
0136 if lukeparam.transpfaste.Dr0 == 0 && lukeparam.transpfaste.Vr0 == 0% cas sans transport, structure vide
0137     transpfaste = '';
0138     transpfaste_id = '';
0139 else
0140     transpfaste = lukeparam.transpfaste;
0141     transpfaste_id = ['_',transpfaste.id];
0142 end
0143 %
0144 % LUKE ripple structure
0145 %
0146 ripple = '';
0147 %
0148 if lukeparam.luke_mode == 2,% save input then exit
0149     %
0150     if auto_mode == 0,
0151         workdir = uigetdir(pwd,'Select save directory');
0152         if ~ischar(workdir)
0153             workdir = pwd;
0154         end
0155     else 
0156         workdir = pwd;
0157     end
0158     %
0159     save([workdir,'/before_luke@',equil_id,'.mat']);
0160     %
0161     return
0162     %
0163 end
0164 %
0165 % LUKE path structure
0166 %
0167 dkepath = getappdata(0,'DKEPATH');
0168 dkepath.temppath = [dkepath.temppathroot,filesep,'CRONOS',filesep,timeid_jd(clock)];%Temporary subdirectory in scratch folder
0169 %
0170 % RUN LUKE with RF waves
0171 %
0172 simul_id = ['METIS_',equil.id,simul_wid,ohm.id,transpfaste_id];
0173 %
0174 dkeparam.psin_S = [];
0175 %
0176 lukestructs = struct;
0177 lukestructs.id_simul = simul_id;
0178 lukestructs.dkeparam = dkeparam;
0179 lukestructs.dkepath = dkepath;
0180 lukestructs.dkedisplay = dkedisplay;
0181 lukestructs.equil = equil;
0182 lukestructs.waves = [];
0183 lukestructs.wavestructs = wavestructs;
0184 lukestructs.ohm = ohm;
0185 lukestructs.transpfaste = transpfaste;
0186 lukestructs.ripple = ripple;
0187 lukestructs.Zf0_interp = [];
0188 %
0189 lukestructs = run_lukert(lukestructs,dkepath);
0190 %
0191 if isfield(lukestructs,'err'),
0192     disp('LUKE run failed for the following reason :')
0193     disperr_jd(lukestructs.err);
0194     error('LUKE run failed')
0195 else
0196     %
0197     Zcurr_RF = lukestructs.output.Zcurr;
0198     ZP0_RF = lukestructs.output.ZP0;
0199     dke_out = lukestructs.output.dke_out;
0200     radialDKE = lukestructs.output.radialDKE;
0201     equilDKE = lukestructs.output.equilDKE;
0202     %
0203 end
0204 %
0205 clear raytracing_yp alphaphimex_jd trapz_dke_yp dmumpsmex
0206 %
0207 signe_out = option.signe;
0208 %
0209 j_TOT_c = signe_out*Zcurr_RF.x_0_vcfsav*mksa.j_ref;
0210 I_TOT_c = sum(signe_out*Zcurr_RF.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref);
0211 p_TOT_c = (ZP0_RF.x_rf_fsav.' + ZP0_RF.x_e_fsav)*mksa.P_ref;
0212 p_RF_c = ZP0_RF.x_rf_fsav.'*mksa.P_ref;
0213 P_TOT_c = sum(p_TOT_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0214 P_RF_c = sum(p_RF_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0215 %
0216 % RUN LUKE without RF waves
0217 %
0218 if sum(abs(ohm.epsi)) ~= 0,%Calculate the Ohmic electric field contribution to the plasma current
0219     %
0220     simul_id = ['METIS_',equil.id,ohm.id,transpfaste_id];
0221     %
0222     lukestructs = struct;
0223     lukestructs.id_simul = simul_id;
0224     lukestructs.dkeparam = dkeparam;
0225     lukestructs.dkepath = dkepath;
0226     lukestructs.dkedisplay = dkedisplay;
0227     lukestructs.equil = equil;
0228     lukestructs.waves = [];
0229     lukestructs.wavestructs = [];
0230     lukestructs.ohm = ohm;
0231     lukestructs.transpfaste = transpfaste;
0232     lukestructs.ripple = ripple;
0233     lukestructs.Zf0_interp = [];
0234     %
0235     lukestructs = run_lukert(lukestructs,dkepath);
0236     %
0237     if isfield(lukestructs,'err'),
0238         disp('LUKE run failed for the following reason :')
0239         disperr_jd(lukestructs.err);
0240         error('LUKE run failed')
0241     else
0242         %
0243         Zcurr = lukestructs.output.Zcurr;
0244         ZP0 = lukestructs.output.ZP0;
0245         %
0246     end
0247     %
0248     j_OHM_c = signe_out*Zcurr.x_0_vcfsav*mksa.j_ref;
0249     I_OHM_c = sum(signe_out*Zcurr.x_0_fsav.*equilDKE.xdA_dke*mksa.j_ref);
0250     p_OHM_c = ZP0.x_e_fsav*mksa.P_ref;
0251     P_OHM_c = sum(p_OHM_c.*equilDKE.xdV_2piRp_dke*2*pi*equilDKE.Rp);
0252 else
0253     j_OHM_c = zeros(1,length(j_TOT_c));
0254     I_OHM_c = 0;
0255     p_OHM_c = zeros(1,length(p_TOT_c)); 
0256     P_OHM_c = 0;
0257 end
0258 %
0259 %j_RF_c = j_TOT_c - j_OHM_c;
0260 I_RF_c = I_TOT_c - I_OHM_c;
0261 %p_OHMsyn_c = p_TOT_c - p_OHM_c - p_RF_c;% synergistic power spent on electric field, already accounted for in E.J
0262 P_OHMsyn_c = P_TOT_c - P_OHM_c - P_RF_c;% synergistic power spent on electric field, already accounted for in E.J
0263 %
0264 disp('-----------------');
0265 disp('   LUKE OUTPUT   ');
0266 disp('-----------------');
0267 disp(' ');
0268 disp(['RF + OHM current  : I_TOT = ',num2str(I_TOT_c,3), 'MA']) 
0269 disp(['OHM current       : I_OHM = ',num2str(I_OHM_c,3), 'MA']) 
0270 disp(['RF current (+syn) : I_RF = ',num2str(I_RF_c,3), 'MA']) 
0271 disp(' ');
0272 disp(['RF + OHM power    : P_TOT = ',num2str(P_TOT_c,3), 'MW']) 
0273 disp(['RF power          : P_RF  = ',num2str(P_RF_c,3), 'MW']) 
0274 disp(['OHM power         : P_OHM = ',num2str(P_OHM_c,3), 'MW']) 
0275 disp(['OHM power (+syn)  : P_OHM = ',num2str(P_OHM_c + P_OHMsyn_c,3), 'MW']) 
0276 %
0277 % interpolation onto metis grid
0278 %
0279 psin = (profil0d.psi - profil0d.psi(1))./(profil0d.psi(end) - profil0d.psi(1));
0280 rho_S = [0,(profil0d.xli(1:end-1) + profil0d.xli(2:end))/2,1];
0281 psin_S = interp1(profil0d.xli,psin,rho_S,'linear');
0282 %
0283 psin_S(1) = 0;% to avoid roundoff errors
0284 psin_S(end) = 1;% to avoid roundoff errors
0285 %
0286 psin_S_mix = unique([psin_S,radialDKE.xpsin_S_dke]);
0287 equilDKE_mix = equilibrium_jd(equil,psin_S_mix);
0288 %
0289 nx = length(equilDKE_mix.xpsin);
0290 nr = length(equilDKE.xpsin);
0291 nm = length(psin);
0292 %
0293 pxpsin_S = ones(nx,1)*radialDKE.xpsin_S_dke;
0294 pxpsin_mix = equilDKE_mix.xpsin.'*ones(1,nr);
0295 %
0296 pxprojmat = (pxpsin_mix > pxpsin_S(:,1:nr)) & (pxpsin_mix < pxpsin_S(:,2:nr+1));
0297 if any(sum(pxprojmat.') ~= ones(1,nx)),
0298     warning('Problem in result interpolation grids')
0299     keyboard
0300 end
0301 %
0302 pppsin_S = ones(nx,1)*psin_S;
0303 pppsin_mix = equilDKE_mix.xpsin.'*ones(1,nm);
0304 %
0305 ppprojmat = (pppsin_mix > pppsin_S(:,1:nm)) & (pppsin_mix < pppsin_S(:,2:nm+1));
0306 if any(sum(ppprojmat.') ~= ones(1,nx)),
0307     warning('Problem in result interpolation grids')
0308     keyboard
0309 end
0310 %
0311 j_TOT_mix = pxprojmat*j_TOT_c.';
0312 j_OHM_mix = pxprojmat*j_OHM_c.';
0313 %
0314 p_TOT_mix = pxprojmat*p_TOT_c.';
0315 p_RF_mix = pxprojmat*p_RF_c.';
0316 p_OHM_mix = pxprojmat*p_OHM_c.';
0317 %
0318 pj_TOT = (ppprojmat.'*(j_TOT_mix.*equilDKE_mix.xdA_dke.'))./(ppprojmat.'*equilDKE_mix.xdA_dke.');
0319 pj_OHM = (ppprojmat.'*(j_OHM_mix.*equilDKE_mix.xdA_dke.'))./(ppprojmat.'*equilDKE_mix.xdA_dke.');
0320 %
0321 pp_TOT = (ppprojmat.'*(p_TOT_mix.*equilDKE_mix.xdV_2piRp_dke.'))./(ppprojmat.'*equilDKE_mix.xdV_2piRp_dke.');
0322 pp_RF = (ppprojmat.'*(p_RF_mix.*equilDKE_mix.xdV_2piRp_dke.'))./(ppprojmat.'*equilDKE_mix.xdV_2piRp_dke.');
0323 pp_OHM = (ppprojmat.'*(p_OHM_mix.*equilDKE_mix.xdV_2piRp_dke.'))./(ppprojmat.'*equilDKE_mix.xdV_2piRp_dke.');
0324 %
0325 output.j_TOT = pj_TOT.'*1e6;
0326 output.j_OHM = pj_OHM.'*1e6;
0327 %
0328 output.I_TOT = I_TOT_c*1e6;
0329 output.I_OHM = I_OHM_c*1e6;
0330 %
0331 output.p_TOT = pp_TOT.'*1e6;
0332 output.p_RF = pp_RF.'*1e6;
0333 output.p_OHM = pp_OHM.'*1e6;
0334 %
0335 output.P_TOT = P_TOT_c*1e6;
0336 output.P_RF = P_RF_c*1e6;
0337 output.P_OHM = P_OHM_c*1e6;
0338 %
0339 output.waves = dke_out.waves;
0340 %
0341 if lukeparam.luke_mode == 1,% save output
0342     %
0343     if auto_mode == 0,
0344         workdir = uigetdir(pwd,'Select save directory');
0345         if ~ischar(workdir)
0346             workdir = pwd;
0347         end
0348     else 
0349         workdir = pwd;
0350     end
0351     %
0352     save([workdir,'/last_luke@',equil_id,'.mat']);
0353     %
0354 end
0355 %
0356 
0357