make_wave_JETliketest

 script make_wave_test_RT

 Parameters for test mode ray-tracing calculations
 This function has a benchmarking purpose only

 by Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr> and J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr>

0001 function [] = make_wave_JETliketest
0002 % script make_wave_test_RT
0003 %
0004 % Parameters for test mode ray-tracing calculations
0005 % This function has a benchmarking purpose only
0006 %
0007 % by Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr> and J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr>
0008 %
0009 close all
0010 %
0011 id_wave = 'JETliketest';
0012 flag_analytic = 2;
0013 %
0014 % Path parameters
0015 %
0016 id_dkepath = '';%For all paths used by DKE solver
0017 path_dkepath = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0018 %
0019 % Equilibrium parameters
0020 %
0021 id_equil = 'JETliketest';%For plasma equilibrium
0022 path_equil = '../EQUIL/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0023 %
0024 % Load structures
0025 %
0026 [equil,dkepath] = load_structures_yp('equil',id_equil,path_equil,'dkepath',id_dkepath,path_dkepath);
0027 %
0028 % =========================================================================
0029 %
0030 % initial ray conditions
0031 %
0032 omega_rf = [3.7]*2*pi*1e9;
0033 %
0034 rho0 = 0.968;
0035 theta0 = 0;
0036 phi0 = 0;
0037 %
0038 m0 = 0;
0039 n0 = NaN;
0040 Nphi0 = -2.0;%initial index of refraction
0041 %
0042 dNpar0 = NaN;
0043 P0_2piRp = NaN;
0044 %
0045 % C3PO computing parameters
0046 %
0047 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)
0048 %
0049 % Display parameters
0050 %
0051 C3POdisplay_pchip.ray = 1;
0052 C3POdisplay_pchip.equilibrium = 0;
0053 C3POdisplay_pchip.p_opt = 2;%Printing or saving option of the figures
0054 C3POdisplay_pchip.mdce = 0;%for distributed computing
0055 %
0056 C3POdisplay_spline.ray = 0;
0057 C3POdisplay_spline.equilibrium = 0;
0058 C3POdisplay_spline.p_opt = -1;%Printing or saving option of the figures
0059 C3POdisplay_spline.mdce = 0;%for distributed computing
0060 %
0061 C3POdisplay_analytic.ray = 1;
0062 C3POdisplay_analytic.equilibrium = 0;
0063 C3POdisplay_analytic.p_opt = -1;%Printing or saving option of the figures
0064 C3POdisplay_analytic.mdce = 0;%for distributed computing
0065 %
0066 % Wave parameters
0067 %
0068 waveparam.mmode = -1;%cold plasma mode [1] : (-1) m (1) p, p is the slow mode when kperp > 0 (ex : LH slow wave)
0069 waveparam.kmode = 0;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0070 %
0071 %Option parameter for FLR effects and cross-comparison between old FP code:
0072 %    - (0): all FLR effects
0073 %    - (1): small FLR effects and 1/vpar dependence
0074 %    - (2): small FLR effects and no 1/vpar dependence and old grid technique for DQL calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0075 %
0076 waveparam.opt_rf = NaN;
0077 %
0078 waveparam.dsmin = NaN;%minimum size for ray fragments
0079 %
0080 % -------------------------------------------------------------------------
0081 %
0082 % Global parameters for the vectorial magnetic equilibrium
0083 %
0084 fitparam.mode_equil = 1;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y)
0085 fitparam.nharm = NaN;%Number of harmonics in the magnetic equilibrium interpolation (NaN, Inf or empty, nharm = ntheta-1)
0086 fitparam.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters
0087 %
0088 % Global parameters for the ray-tracing
0089 %
0090 rayparam.testmode = 1;
0091 rayparam.tensortype = waveparam.kmode;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0092 rayparam.t0 = 0;
0093 rayparam.tfinal = 30000;
0094 rayparam.dt0 = 1.e-4;
0095 rayparam.dS = 1.e-4;
0096 rayparam.tol = 1e-12;%when tolerance is increased (less accurate calculation of D=0), tfinal must be increased accordingly
0097 rayparam.kmax = 30000;
0098 rayparam.ncyclharm = 3;%number of cyclotron harmonics (just for hot and relativistic dielectric tensors)
0099 rayparam.reflection = 1;%1:Enforce wave reflection at plasma boundary, 0: the code calculates itself if the ray must leave of not the plasma
0100 rayparam.rel_opt = 1;%option for (1) relativistic or (0) non-relativistic calculations
0101 rayparam.nperp = 1000;%number of points in pperp integration for damping calculations
0102 rayparam.pperpmax = 10;%maximum value of pperp in damping calculations
0103 rayparam.tau_lim = 20;%value of optical depth beyond which the wave is considered absorbed
0104 rayparam.metricmode = 1;% 1 : Only curvilinear metric; 2 : Both curvilinear and cartesian metric
0105 rayparam.rhoswitch = 0.1;% Radial position at which the switch occurs when metricmode = 2
0106 rayparam.substitution_method = 7;%Method of substitution : 1->dichotomy1D, 2->dichotomy2D, 3->newton_numeric1D, 4->newton_numeric2D, 5->newton_analytic1D, 6->newton_analytic2D, 7->newton_mix2D
0107 %
0108 % =========================================================================
0109 %
0110 % C3P0 ray tracing
0111 %
0112 fitparam.method = 'pchip';%nearest,spline,pchip
0113 tstart = tic;
0114 equil_fit_pchip = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);%Build vectorized magnetic equilibrium structure
0115 telapsed_equil_pchip = toc(tstart);
0116 %
0117 info_dke_yp(2,['Vectorial form of the magnetic equilibrium ',equil.id,' is calculated with pchip method.']);
0118 if C3POdisplay_pchip.equilibrium,testfitequil_yp(equil,equil_fit_pchip);end
0119 %
0120 fitparam.method = 'spline';%nearest,spline,pchip
0121 tstart = tic;
0122 equil_fit_spline = fitequil_yp(equil,fitparam.mode_equil,fitparam.method,fitparam.ngridresample,fitparam.nharm);%Build vectorized magnetic equilibrium structure
0123 telapsed_equil_spline = toc(tstart);
0124 %
0125 info_dke_yp(2,['Vectorial form of the magnetic equilibrium ',equil.id,' is calculated with spline method.']);
0126 if C3POdisplay_spline.equilibrium,testfitequil_yp(equil,equil_fit_spline);end
0127 %
0128 % Calculation of Npar0 using Nphi0
0129 %
0130 % warning : this is valid only for m=0 cases
0131 %
0132 Bx_a0_fit = ppval(equil_fit_pchip.Bx_fit.pp_a0,rho0);
0133 Bx_an_fit = ppval(equil_fit_pchip.Bx_fit.pp_an,rho0);
0134 Bx_bn_fit = ppval(equil_fit_pchip.Bx_fit.pp_bn,rho0);
0135 %
0136 By_a0_fit = ppval(equil_fit_pchip.By_fit.pp_a0,rho0);
0137 By_an_fit = ppval(equil_fit_pchip.By_fit.pp_an,rho0);
0138 By_bn_fit = ppval(equil_fit_pchip.By_fit.pp_bn,rho0);
0139 %
0140 Bz_a0_fit = ppval(equil_fit_pchip.Bz_fit.pp_a0,rho0);
0141 Bz_an_fit = ppval(equil_fit_pchip.Bz_fit.pp_an,rho0);
0142 Bz_bn_fit = ppval(equil_fit_pchip.Bz_fit.pp_bn,rho0);
0143 %
0144 B_a0_fit = ppval(equil_fit_pchip.B_fit.pp_a0,rho0);
0145 B_an_fit = ppval(equil_fit_pchip.B_fit.pp_an,rho0);
0146 B_bn_fit = ppval(equil_fit_pchip.B_fit.pp_bn,rho0);
0147 %
0148 % Build interpolated magnetic fields
0149 %
0150 [xBx] = calcval_yp(equil_fit_pchip,theta0,Bx_a0_fit,Bx_an_fit,Bx_bn_fit);
0151 [xBy] = calcval_yp(equil_fit_pchip,theta0,By_a0_fit,By_an_fit,By_bn_fit);
0152 [xBz] = calcval_yp(equil_fit_pchip,theta0,Bz_a0_fit,Bz_an_fit,Bz_bn_fit);
0153 [xB] = calcval_yp(equil_fit_pchip,theta0,B_a0_fit,B_an_fit,B_bn_fit);
0154 %
0155 xBzn = xBz./xB;
0156 %
0157 Npar0 = Nphi0.*xBzn;
0158 %
0159 rayinit.omega_rf = omega_rf;
0160 rayinit.yrho0 = rho0;%Initial radial position at launch
0161 rayinit.ytheta0 = theta0;%Initial poloidal position at launch
0162 rayinit.yphi0 = phi0;%Initial axial (toroidal) position at launch
0163 rayinit.ym0 = m0;%Initial poloidal mode number
0164 rayinit.yn0 = n0;%Initial wave vector along the axial (toroidal) direction
0165 rayinit.yNpar0 = Npar0;%Initial index of refraction
0166 rayinit.ydNpar0 = dNpar0;%initial Ray spectral width
0167 rayinit.yP0_2piRp = P0_2piRp;%Lineic initial power density initial power in the ray (W/m)
0168 %
0169 % --------------------------------------------------------------------------
0170 %
0171 % C3PO computing parameters
0172 %
0173 C3POparam.clustermode.main_C3PO_jd.scheduler.mode = mdce_mode_main_C3PO_jd;%MatLab distributed computing environment
0174 %
0175 % ray-tracing calculations
0176 %
0177 tstart = tic;
0178 wave_numeric_pchip = main_C3PO_jd(dkepath,[id_wave,'_numeric_pchip'],equil,equil_fit_pchip,rayinit,waveparam,[],rayparam,C3POdisplay_pchip,C3POparam,[],[],0);clear mex;clear functions
0179 telapsed_ray_numeric_pchip = toc(tstart);
0180 %
0181 %info_dke_yp(2,'Ray trajectories calculated (interpolated magnetic equilibrium with pchip method)');
0182 %
0183 % Calculation of Npar0 using Nphi0
0184 %
0185 % warning : this is valid only for m=0 cases
0186 %
0187 Bx_a0_fit = ppval(equil_fit_spline.Bx_fit.pp_a0,rho0);
0188 Bx_an_fit = ppval(equil_fit_spline.Bx_fit.pp_an,rho0);
0189 Bx_bn_fit = ppval(equil_fit_spline.Bx_fit.pp_bn,rho0);
0190 %
0191 By_a0_fit = ppval(equil_fit_spline.By_fit.pp_a0,rho0);
0192 By_an_fit = ppval(equil_fit_spline.By_fit.pp_an,rho0);
0193 By_bn_fit = ppval(equil_fit_spline.By_fit.pp_bn,rho0);
0194 %
0195 Bz_a0_fit = ppval(equil_fit_spline.Bz_fit.pp_a0,rho0);
0196 Bz_an_fit = ppval(equil_fit_spline.Bz_fit.pp_an,rho0);
0197 Bz_bn_fit = ppval(equil_fit_spline.Bz_fit.pp_bn,rho0);
0198 %
0199 B_a0_fit = ppval(equil_fit_spline.B_fit.pp_a0,rho0);
0200 B_an_fit = ppval(equil_fit_spline.B_fit.pp_an,rho0);
0201 B_bn_fit = ppval(equil_fit_spline.B_fit.pp_bn,rho0);
0202 %
0203 % Build interpolated magnetic fields
0204 %
0205 [xBx] = calcval_yp(equil_fit_spline,theta0,Bx_a0_fit,Bx_an_fit,Bx_bn_fit);
0206 [xBy] = calcval_yp(equil_fit_spline,theta0,By_a0_fit,By_an_fit,By_bn_fit);
0207 [xBz] = calcval_yp(equil_fit_spline,theta0,Bz_a0_fit,Bz_an_fit,Bz_bn_fit);
0208 [xB] = calcval_yp(equil_fit_spline,theta0,B_a0_fit,B_an_fit,B_bn_fit);
0209 %
0210 xBzn = xBz./xB;
0211 %
0212 Npar0 = Nphi0.*xBzn;
0213 %
0214 rayinit.omega_rf = omega_rf;
0215 rayinit.yrho0 = rho0;%Initial radial position at launch
0216 rayinit.ytheta0 = theta0;%Initial poloidal position at launch
0217 rayinit.yphi0 = phi0;%Initial axial (toroidal) position at launch
0218 rayinit.ym0 = m0;%Initial poloidal mode number
0219 rayinit.yn0 = n0;%Initial wave vector along the axial (toroidal) direction
0220 rayinit.yNpar0 = Npar0;%Initial index of refraction
0221 rayinit.ydNpar0 = dNpar0;%initial Ray spectral width
0222 rayinit.yP0_2piRp = P0_2piRp;%Lineic initial power density initial power in the ray (W/m)
0223 %
0224 % C3PO computing parameters
0225 %
0226 C3POparam.clustermode.main_C3PO_jd.scheduler.mode = mdce_mode_main_C3PO_jd;%MatLab distributed computing environment
0227 %
0228 % ray-tracing calculations
0229 %
0230 %tstart = tic;
0231 %wave_numeric_spline = main_C3PO_jd(dkepath,[id_wave,'_numeric_spline'],equil,equil_fit_spline,rayinit,waveparam,[],rayparam,C3POdisplay_spline,C3POparam,[],[],0);clear mex;clear functions
0232 %telapsed_ray_numeric_spline = toc(tstart);
0233 %
0234 info_dke_yp(2,'Ray trajectories calculated (interpolated magnetic equilibrium with spline method)');
0235 %
0236 %tstart = tic;
0237 %wave_analytic = main_C3PO_jd(dkepath,[id_wave,'_analytic'],equil,equil_fit_pchip,rayinit,waveparam,[],rayparam,C3POdisplay_analytic,C3POparam,[],[],flag_analytic);
0238 %telapsed_ray_analytic = toc(tstart);
0239 %
0240 %info_dke_yp(2,'Ray trajectories calculated (analytic magnetic equilibrium)');
0241 %
0242 save_str = ['WAVE_',id_wave,'.mat'];
0243 save(save_str,'id_wave','wave_numeric_pchip','equil_fit_pchip');
0244 %
0245 %info_dke_yp(2,'Wave parameters saved');
0246 %
0247 disp(['Elapsed time (s):',num2str(telapsed_equil_pchip),' for interpolating MHD equilibrium (PCHIP interp. technique)']);
0248 %disp(['Elapsed time (s):',num2str(telapsed_equil_spline),' for interpolating MHD equilibrium (SPLINE interp. technique)']);
0249 %
0250 disp(['Elapsed time (s):',num2str(telapsed_ray_numeric_pchip),' for ray calculations with numerical MHD equilibrium (PCHIP interp. technique)']);
0251 %disp(['Elapsed time (s):',num2str(telapsed_ray_numeric_spline),' for ray calculations with numerical MHD equilibrium (SPLINE interp. technique)']);
0252 %disp(['Elapsed time (s):',num2str(telapsed_ray_analytic),' for ray calculations with analytical MHD equilibrium']);
0253 %
0254 % --- display results ---
0255 %
0256 waves = {wave_numeric_pchip};
0257 rays = {wave_numeric_pchip.rays{1}};
0258 %
0259 legs = {'Numeric - pchip';...
0260        'Numeric_pchip'};
0261 %
0262 equil_fit = [equil_fit_pchip];
0263 %
0264 filename = ['Fig_',id_wave];
0265 %
0266 opt.p_opt = C3POdisplay_pchip.p_opt;
0267 opt.ntheta_fit = 65;
0268 opt.nrho_fit = 15;
0269 opt.propvar = 1;   
0270 %
0271 graph_comp_RT_jd(rays,legs,equil_fit,filename,opt)
0272 %
0273 diary4cvs_C3PO_yp(id_wave,dkepath,waves);% diary some results for CVS validation