make_wave_JETliketest_cyl

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