make_wave_JETliketest_bis

 script make_wave_test_RT

 Parameters for test mode ray-tracing calculations
 This function is intended fr benchmarking purpose or regression tests

   Inputs :
   
      - mode: purpose of the script, (0) full benchmark, (>0) regression tests [1,1]

   Outputs : 

      - waves: wave structures
      - rays: rays structures
      - equils_fit: fitted equilibria structures
     
 by Y. Peysson (DRFC/DSM/CEA) <yves.peysson@cea.fr> and J. Decker (DRFC/DSM/CEA) <joan.decker@cea.fr>

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