rundke_ITER_Scen2_C3PO_EC_fluctn_sigma0p1_delta0p01

 Script for running LUKE

 This script is automatically generated by irunluke_jd.
 by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)

0001 % Script for running LUKE
0002 %
0003 % This script is automatically generated by irunluke_jd.
0004 % by J. Decker (joan.decker@cea.fr) and Y. Peysson (yves.peysson@cea.fr)
0005 %
0006 clear all
0007 clear mex
0008 clear functions
0009 close all
0010 warning('off','all')
0011 pause(1)
0012 %
0013 p_opt = -1;
0014 %
0015 % *********************** Specify LUKE structures *****************************
0016 %
0017 id_simul = '';%Simulation ID
0018 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0019 %
0020 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0021 rho_S = [0:0.1:0.6,0.61:0.01:0.79,0.80:0.05:1];
0022 %
0023 id_dkepath = '';%For all paths used by DKE solver
0024 path_dkepath = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0025 %
0026 id_equil = 'ITER_Scen2_200103121816_430_129';%For plasma equilibrium
0027 path_equil = 'EQUIL/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0028 %
0029 id_fluct = 'ne_sigma0p1_delta0p01';%For fast electron radial transport
0030 path_fluct = 'FLUCT/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0031 %
0032 id_dkeparam = 'EC_RT';%For DKE code parameters
0033 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_dkedisplay = 'NO_DISPLAY';%For output code display
0036 path_dkedisplay = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_ohm = '';%For Ohmic electric contribution
0039 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0042 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 id_transpfaste = '';%For fast electron radial transport
0045 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_ripple = '';%For fast electron magnetic ripple losses
0048 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 %************************************************************************************************************************************
0051 %
0052 % LAUNCH parameters
0053 %
0054 id = 'NTM_32_fluctn';
0055 freq_GHz = 170;
0056 mmode = -1;% cold plasma mode (1) m (-1) p, (slow/fast LH wave: -1/+1, O/X mode: -1/+1) :1
0057 %
0058 R_L = 6.4848;
0059 Z_L = 4.110;
0060 P_L = 20;% [MW]
0061 %
0062 cone_div = 1.08;% equivalent exp(-2) power cone divergence
0063 cone_ori = -0.5;% equivalent cone origin
0064 %
0065 z_L = 0.35;% focal point
0066 %
0067 alpha_prater = -63.5;% angle between N_P and horizontal plane
0068 beta_prater = 22;% angle between N and poloidal plane
0069 %
0070 ns = 1000;
0071 method = 'cubic';
0072 %
0073 % C3PO computing parameters
0074 %
0075 mdce_mode_main_C3PO_jd = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2)for the function main_C3PO_jd.m (MDC toolbox must be installed for option 1)
0076 %
0077 % Display parameters
0078 %
0079 C3POdisplay.ray = 0;
0080 C3POdisplay.equilibrium = 0;
0081 C3POdisplay.p_opt = -1;%Printing or saving option of the figures
0082 %
0083 % Wave parameters
0084 %
0085 waveparam.mmode = mmode;%cold plasma mode [1] : (-1) m (1) p, p is the slow mode when kperp > 0 (ex : LH slow wave)
0086 waveparam.kmode = 0;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0087 %
0088 %Option parameter for FLR effects and cross-comparison between old FP code:
0089 %    - (0): all FLR effects
0090 %    - (1): small FLR effects and 1/vpar dependence
0091 %    - (2): small FLR effects and no 1/vpar dependence and old grid technique for DQL calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0092 %
0093 waveparam.opt_rf = 1;
0094 %
0095 waveparam.dsmin = 0;%minimum size for ray fragments
0096 %
0097 waveparam.nd = 1;%Number of transverse distance positions for beamlets
0098 waveparam.nchi = 1;%Number of angular positions for beamlets
0099 %
0100 waveparam.n_rf_list = 1:2;
0101 waveparam.ns = 1;%ray smoothing parameter
0102 %
0103 % -------------------------------------------------------------------------
0104 %
0105 % Global parameters for the vectorial magnetic equilibrium
0106 %
0107 % Equilibrium parameters for the ray-tracing
0108 %
0109 fitparam.mode_equil = 1;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y)
0110 fitparam.method = 'spline';%nearest,spline,pchip,cubic
0111 fitparam.nharm = NaN;%Number of harmonics in the magnetic equilibrium interpolation (less than ntheta_equil/2)
0112 fitparam.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters
0113 fitparam.opt_load = 0;%Reload existing vectorial magnetic equilibrium (1) or overwrite it (0).
0114 %
0115 % Global parameters for the ray-tracing
0116 %
0117 rayparam.testmode = 0;
0118 rayparam.tensortype = 0;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0119 rayparam.t0 = 0;
0120 rayparam.tfinal = 50000;
0121 rayparam.dt0 = 0.0001;
0122 rayparam.dS = 0.001;
0123 rayparam.tol = 1e-8;%when tolerance is increased (less accurate calculation of D=0), tfinal must be increased accordingly
0124 rayparam.kmax = 50000;
0125 rayparam.ncyclharm = 3;%number of cyclotron harmonics (just for hot and relativistic dielectric tensors)
0126 rayparam.reflection = 0;%1:Enforce wave reflection at plasma boundary, 0: the code calculates itself if the ray must leave of not the plasma
0127 rayparam.rel_opt = 1;%option for (1) relativistic or (0) non-relativistic calculations
0128 rayparam.nperp =1000;%number of points in pperp integration for damping calculations
0129 rayparam.pperpmax = 10;%maximum value of pperp in damping calculations
0130 rayparam.tau_lim = 20;%value of optical depth beyond which the wave is considered absorbed
0131 rayparam.kextra = 1000;
0132 %
0133 %************************************************************************************************************************************
0134 %
0135 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,fluct] = load_structures_yp('dkepath',id_dkepath,path_dkepath,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_dkedisplay,path_dkedisplay,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple,'fluct',id_fluct,path_fluct);
0136 %
0137 % Parallel computing parameters
0138 %
0139 C3POparam.clustermode.main_C3PO_jd.scheduler.mode = mdce_mode_main_C3PO_jd;%MatLab distributed computing environment
0140 %
0141 % launch structure
0142 %
0143 launch.id = id;
0144 launch.type = 'EC';
0145 launch.omega_rf = [freq_GHz]*2*pi*1e9;%wave angular frequency (rad/s)
0146 %
0147 launch.yR_L = R_L;
0148 launch.yZ_L = Z_L;
0149 launch.yphi_L = 0;
0150 launch.yalpha_L = sign(beta_prater)*pi - atan(tand(beta_prater)/cosd(alpha_prater));
0151 launch.ybeta_L = acos(cosd(beta_prater)*sind(alpha_prater));
0152 launch.yP_L = P_L*1e6;
0153 launch.dNpar0 = NaN;
0154 launch.ns = ns;
0155 launch.method = method;
0156 %
0157 launch.w0 = 0.3/(freq_GHz*pi*cone_div*pi/180);
0158 launch.z_L = z_L;
0159 %
0160 % wavestruct
0161 %
0162 wavestruct.wavesolver = 'C3PO';
0163 wavestruct.launch = launch;
0164 wavestruct.id = [wavestruct.wavesolver,'_',launch.id];
0165 wavestruct.waveparam = waveparam;
0166 wavestruct.rayparam = rayparam;
0167 wavestruct.fitequilparam = fitparam;
0168 wavestruct.raydisplay = C3POdisplay;
0169 wavestruct.C3POparam = C3POparam;
0170 %
0171 id_simul = [id_equil,'_',wavestruct.id,'_',id_fluct];
0172 %
0173 fluct_fitparam.method = 'pchip';%nearest,spline,pchip,cubic
0174 fluct_fitparam.nharm = 32;%Number of harmonics in the plasma fluctuations interpolation (less than ntheta_equil/2)
0175 fluct_fitparam.ngridresample = 201;%Number of grid points for resampling the radial profile of plasma fluctuations parameters (very slow if too big !!)
0176 %
0177 %************************************************************************************************************************************
0178 %
0179 dkeparam.rt_mode = 1;
0180 dkeparam.dtn = 1i;% as an imaginary number, dkeparam.dtn prescribes the number of internal time steps
0181 dkeparam.tn = 20i;% as an imaginary number, dkeparam.tn prescribes the number of fluctuation iterations
0182 dkeparam.Nfluct = 1.052872583699049;% to retrieve previous results before merge
0183 %
0184 dkedisplay.display_mode = 1;
0185 dkedisplay.display_time_mode = 1;
0186 %
0187 opt_save = 1;%In case of crash, to restart simulation from an intermediate step stored in a *.mat file (<0 for launching simulation from backup, 0 no backup, >0 for performing backup)
0188 %
0189 equil.fluct = fluct;
0190 equil.fluct.fitparam = fluct_fitparam;
0191 %
0192 if exist('dmumpsmex','file');
0193     dkeparam.invproc = -2;%MUMPSMEX
0194 end
0195 %
0196 %************************************************************************************************************************************
0197 %
0198 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = run_lukert(id_simul,{wavestruct},[],psin_S,rho_S,dkepath,equil,dkeparam,dkedisplay,ohm,transpfaste,ripple,[],[],opt_save);
0199 %
0200 % ---------------------------------------------- Save results ----------------------------------------------
0201 %
0202 filename = [path_simul,'LUKE_RESULTS_',id_simul,'.mat'];
0203 save(filename);
0204 info_dke_yp(2,['LUKE Results saved in ',filename]);
0205 %
0206 % ---------------------------------------------- Display results ----------------------------------------------
0207 %
0208 p_opt = 2;
0209 path_simul = pwd;
0210 [tn,tn_all,fluctres] = fluctanalysis_yp(equilDKE,ZP0,Zcurr,dkeparam,mksa,0);
0211 %
0212 info_dke_yp(2,['Total co-current     : ',num2str(fluctres.I.Itotm*fluctres.I.scocurr),' +/- ',num2str(abs(fluctres.I.Itots)*fluctres.I.scocurr),' MA, relative unbiased accuracy: ',num2str(abs(fluctres.I.Itotr))]);
0213 info_dke_yp(2,['Total RF power    : ',num2str(fluctres.P.Ptotm),' +/- ',num2str(abs(fluctres.P.Ptots)),' MW, relative unbiased accuracy: ',num2str(abs(fluctres.P.Ptotr))]);
0214 disp(['Rejected simulations (Pin not equal to Pabs within error bars): ',int2str(fluctres.P.iPtotr(:)')]);
0215 %
0216 figure('Name','Time evolution EC current (fluctuation)')
0217 %
0218 if ~isscalar(fluctres.I.xJs), 
0219     graph1D_jd(fluctres.I.tnI,fluctres.I.tI*fluctres.I.scocurr,0,0,'t\nu_{coll}','I_{lh} co-current (MA)',[],NaN,'','','-','none','r',2,20,gca,1,NaN);
0220     graph1D_jd(fluctres.I.tnI,fluctres.I.tI*fluctres.I.scocurr,0,0,'t\nu_{coll}','I_{lh} co-current (MA)',[],NaN,'','','-','none','r',2,20,gca,1,NaN);
0221     graph1D_jd(fluctres.I.tnI(fluctres.I.iItot),fluctres.I.tI(fluctres.I.iItot)*fluctres.I.scocurr,0,0,'t\nu_{coll}','I_{lh} co-current (MA)',[],NaN,'','','none','+','r',2,20,gca,1,NaN);
0222     graph1D_jd(fluctres.I.tnI(fluctres.I.iItot),fluctres.I.Itotm*fluctres.I.scocurr*ones(1,length(tn(fluctres.I.iItot))),0,0,'t\nu_{coll}','I_{lh} (MA)',['I_{tot} = ',num2str(fluctres.I.Itotm),'\pm',num2str(fluctres.I.Itots*fluctres.I.scocurr),' (MA)'],NaN,'','','--','none','r',1,20,gca,1,NaN);
0223     graph1D_jd(fluctres.I.tnI(fluctres.I.iItot),fluctres.I.Itotm*fluctres.I.scocurr*ones(1,length(fluctres.I.iItot)),0,0,'t\nu_{coll}','I_{lh} co-current (MA)',['I_{tot} = ',num2str(fluctres.I.Itotm),'\pm',num2str(fluctres.I.Itots*fluctres.I.scocurr),' (MA)'],NaN,[min(tn),max(tn)],'','none','+','r',1,20,gca,1,NaN,NaN,fluctres.I.Itots*ones(1,length(fluctres.I.iItot)),fluctres.I.Itots*ones(1,length(fluctres.I.iItot)));
0224     %legend(['t_{fluct}\nu_{coll} = ',num2str(dtn_fluct)],['HWHM_{fluct} = ',num2str(hwhm_npar)]);
0225     %
0226     print_jd(p_opt,[id_simul,'_IEC_time'],path_simul);
0227 end
0228 %
0229 figure('Name','Time evolution EC power (fluctuation)')
0230 %
0231 if ~isscalar(fluctres.P.xPs), 
0232     graph1D_jd(tn,fluctres.P.tP,0,0,'t*\nu_{coll}','P_{lh} (MW)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,[min(tn),max(tn)],[fluctres.P.Ptotm-10*abs(fluctres.P.Ptots),fluctres.P.Ptotm+10*abs(fluctres.P.Ptots)],'-','none','b',2,20,gca,1,NaN);
0233     graph1D_jd(tn(fluctres.P.iPtot),fluctres.P.tP(fluctres.P.iPtot),0,0,'t\nu_{coll}','P_{lh} (MW)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,[min(tn),max(tn)],[fluctres.P.Ptotm-10*abs(fluctres.P.Ptots),fluctres.P.Ptotm+10*abs(fluctres.P.Ptots)],'+','none','b',2,20,gca,1,NaN);
0234     graph1D_jd(tn,fluctres.P.Ptotm*ones(1,length(tn)),0,0,'t\nu_{coll}','P_{lh} (MW)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,[min(tn),max(tn)],'','--','none','b',1,20,gca,1,NaN);
0235     graph1D_jd(tn(fluctres.P.iPtot),fluctres.P.Ptotm*ones(1,length(fluctres.P.iPtot)),0,0,'t\nu_{coll}','P_{lh} (MW)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,[min(tn),max(tn)],'','none','+','b',1,20,gca,1,NaN,NaN,fluctres.P.Ptots*ones(1,length(fluctres.P.iPtot)),fluctres.P.Ptots*ones(1,length(fluctres.P.iPtot)));
0236     %
0237     print_jd(p_opt,[id_simul,'_PEC_time'],path_simul);
0238 end
0239 %
0240 figure('Name','Mean co-current profile (fluctuation)')
0241 %
0242 index = find(fluctres.I.xJm*fluctres.I.scocurr > 1e-5*max(fluctres.I.xJm*fluctres.I.scocurr))
0243 [curve, goodness] = fit(equilDKE.xrho(index)',log(abs(fluctres.I.xJm(index)*fluctres.I.scocurr))','poly2');
0244 fitmax = max(exp(curve.p1*equilDKE.xrho(index).^2 + curve.p2*equilDKE.xrho(index) + curve.p3));
0245 xJmmax = max(fluctres.I.xJm*fluctres.I.scocurr);
0246 rhomax = equilDKE.xrho(fluctres.I.xJm*fluctres.I.scocurr == xJmmax);
0247 delta_Jm = 2*sqrt(-log(2)/curve.p1);
0248 %
0249 if isscalar(fluctres.I.xJs), 
0250     graph1D_jd(equilDKE.xrho,fluctres.I.xJm*fluctres.I.scocurr,0,0,'\rho','J_{lh} co-current (MA/m^2)',['I_{lh} (co-current) = ',num2str(fluctres.I.Itotm*fluctres.I.scocurr),'\pm',num2str(abs(fluctres.I.Itots)*fluctres.I.scocurr),' (MA)'],NaN,'','','-','none','r',2,20,gca,1,NaN);
0251 else
0252     graph1D_jd(equilDKE.xrho,fluctres.I.xJm*fluctres.I.scocurr,0,0,'\rho','J_{lh} co-current (MA/m^2)',['I_{lh} (co-current) = ',num2str(fluctres.I.Itotm*fluctres.I.scocurr),'\pm',num2str(abs(fluctres.I.Itots)*fluctres.I.scocurr),' (MA)'],NaN,'','','-','none','r',1,20,gca,1,NaN,NaN,fluctres.I.xJs*fluctres.I.scocurr,fluctres.I.xJs*fluctres.I.scocurr);
0253     graph1D_jd(equilDKE.xrho,fluctres.I.xJm*fluctres.I.scocurr,0,0,'\rho','J_{lh} co-current (MA/m^2)',['I_{lh} (co-current) = ',num2str(fluctres.I.Itotm*fluctres.I.scocurr),'\pm',num2str(abs(fluctres.I.Itots)*fluctres.I.scocurr),' (MA)'],NaN,'','','-','none','r',2,20,gca,1,NaN);
0254     graph1D_jd(equilDKE.xrho(index),xJmmax*exp(curve.p1*equilDKE.xrho(index).^2 + curve.p2*equilDKE.xrho(index) + curve.p3)/fitmax,0,0,'\rho','J_{lh} co-current (MA/m^2)',['I_{lh} (co-current) = ',num2str(fluctres.I.Itotm*fluctres.I.scocurr),'\pm',num2str(abs(fluctres.I.Itots)*fluctres.I.scocurr),' (MA)'],NaN,'','','--','none','k',1,20,gca,1,NaN);
0255 end
0256 %
0257 disp(['FWHM of the current density profile = ',num2str(delta_Jm)]);
0258 disp(['Position of the current density profile = ',num2str(rhomax)]);
0259 %
0260 print_jd(p_opt,[id_simul,'_jcocurr_EC_profile'],path_simul);
0261 %
0262 figure('Name','Mean EC power absorption profile (fluctuation)')
0263 %
0264 if isscalar(fluctres.P.xPs),
0265     graph1D_jd(equilDKE.xrho,fluctres.P.xPm,0,0,'\rho','P_{tot} (MW/m^3)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,'','','-','none','b',2,20,gca,1,NaN);
0266 else
0267     graph1D_jd(equilDKE.xrho,fluctres.P.xPm,0,0,'\rho','P_{tot} (MW/m^3)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,'','','-','none','b',1,20,gca,1,NaN,NaN,fluctres.P.xPs,fluctres.P.xPs);
0268     graph1D_jd(equilDKE.xrho,fluctres.P.xPm,0,0,'\rho','P_{tot} (MW/m^3)',['P_{tot} = ',num2str(fluctres.P.Ptotm),'\pm',num2str(fluctres.P.Ptots),' (MW)'],NaN,'','','-','none','b',2,20,gca,1,NaN);
0269 end
0270 print_jd(p_opt,[id_simul,'_p_EC_profile'],path_simul);