make_wave_JETlikeflucttest_noripple

 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_noripple
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 = 'JETlikeflucttest_noripple';
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';%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 % Density and magnetic field fluctuation
0024 %
0025 id_fluct = 'noripple';%For density and magnetic field fluctuation
0026 path_fluct = '../FLUCT/';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0027 %
0028 % Load structures
0029 %
0030 [equil,dkepath,fluct] = load_structures_yp('equil',id_equil,path_equil,'dkepath',id_dkepath,path_dkepath,'fluct',id_fluct,path_fluct);
0031 %
0032 % =========================================================================
0033 %
0034 % initial ray conditions
0035 %
0036 omega_rf = [3.7]*2*pi*1e9;
0037 %
0038 rho0 = 0.968;
0039 theta0 = [pi/12];
0040 phi0 = 0;%between two toroidal magnetic field coils (phi = 0 corresponds to a coil position)
0041 %
0042 m0 = 0;
0043 n0 = NaN;
0044 Nphi0 = -2.0;%initial index of refraction
0045 %
0046 dNpar0 = NaN;
0047 P0_2piRp = NaN;
0048 %
0049 % C3PO computing parameters
0050 %
0051 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)
0052 %
0053 % Display parameters
0054 %
0055 C3POdisplay.ray = 0;
0056 C3POdisplay.equilibrium = 0;
0057 C3POdisplay.fluctuations = 0;
0058 C3POdisplay.p_opt = 2;%Printing or saving option of the figures
0059 C3POdisplay.mdce = 1;%for distributed computing
0060 %
0061 % Wave parameters
0062 %
0063 waveparam.mmode = -1;%cold plasma mode [1] : (-1) m (1) p, p is the slow mode when kperp > 0 (ex : LH slow wave)
0064 waveparam.kmode = 0;%(0:cold,1:warm,2:hot;3:weak realtivistic,4:full relativistic)
0065 %
0066 %Option parameter for FLR effects and cross-comparison between old FP code:
0067 %    - (0): all FLR effects
0068 %    - (1): small FLR effects and 1/vpar dependence
0069 %    - (2): small FLR effects and no 1/vpar dependence and old grid technique for DQL calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0070 %
0071 waveparam.opt_rf = NaN;
0072 %
0073 waveparam.dsmin = NaN;%minimum size for ray fragments
0074 %
0075 % -----------------------------------------------------------------------------------------------
0076 %
0077 % Global parameters for the vectorial magnetic equilibrium and the plasma fluctuations (if calculated)
0078 %
0079 fitparam.equil.method = 'spline';%nearest,spline,pchip
0080 fitparam.equil.nharm = NaN;%Number of harmonics in the magnetic equilibrium interpolation (NaN, Inf or empty, nharm = ntheta-1)
0081 fitparam.equil.ngridresample = 1001;%Number of grid points for resampling the radial profile of magnetic equilibrium parameters
0082 fitparam.equil.mode_equil = 1;%(rho,theta) -> 0, (x,y) -> 1
0083 %
0084 fitparam.fluct.mode_equil = fitparam.equil.mode_equil;%Magnetic equilibrium grid type: (1): (psi-theta), (2): (x-y)
0085 fitparam.fluct.method = 'pchip';%nearest,spline,pchip
0086 fitparam.fluct.nharm = 32;%Number of harmonics in the plasma fluctuations interpolation (less than ntheta_equil/2)
0087 fitparam.fluct.ngridresample = 201;%Number of grid points for resampling the radial profile of plasma fluctuations parameters (very slow if too big !!)
0088 %
0089 % Global parameters for the ray-tracing
0090 %
0091 rayparam.testmode = 0;
0092 rayparam.tensortype = waveparam.kmode;%(0:cold,1:warm,2:hot;3:weak relativistic,4:full relativistic)
0093 rayparam.t0 = 0;
0094 rayparam.tfinal = 10000;
0095 rayparam.dt0 = 1.e-4;
0096 rayparam.dS = 1.e-4;
0097 rayparam.tol = 1e-12;%when tolerance is increased (less accurate calculation of D=0), tfinal must be increased accordingly
0098 rayparam.kmax = 60000;%
0099 rayparam.ncyclharm = 3;%number of cyclotron harmonics (just for hot and relativistic dielectric tensors)
0100 rayparam.reflection = 1;%1:Enforce wave reflection at plasma boundary, 0: the code calculates itself if the ray must leave of not the plasma
0101 rayparam.rel_opt = 1;%option for (1) relativistic or (0) non-relativistic calculations
0102 rayparam.nperp = 10000;%number of points in pperp integration for damping calculations
0103 rayparam.pperpmax = 10;%maximum value of pperp in damping calculations
0104 rayparam.tau_lim = 20;%value of optical depth beyond which the wave is considered absorbed (usually 20. Otherwise Inf)
0105 rayparam.kextra = 1000;%number of calculations performed beyond the full linear absorption (for quasilinear calculations which may require more points)
0106 %
0107 % -------------------------------------------------------------------------
0108 %
0109 % C3P0 ray tracing
0110 %
0111 % Vectorial description of the magnetic equilibrium
0112 %
0113 equil_fit = fitequil_yp(equil,fitparam.equil.mode_equil,fitparam.equil.method,fitparam.equil.ngridresample,fitparam.equil.nharm);%Build vectorized magnetic equilibrium structure
0114 info_dke_yp(2,['Vectorial form of the magnetic equilibrium ',equil.id,' is calculated.']);
0115 if C3POdisplay.equilibrium,testfitequil_yp(equil,equil_fit);end
0116 %
0117 % Vectorial description of the plasma fluctuations
0118 %
0119 if ~isempty(fluct),
0120    fluct = fluctphase_yp(fluct,0);%Mersenne-Twitter pseudorandom generator -> 0, Matlab rand (reset by clock) -> 1
0121    [fluct_fit] = fitfluct_yp(fluct,fitparam.equil.mode_equil,fitparam.fluct.method,fitparam.fluct.ngridresample,fitparam.fluct.nharm);%Build vectorized plasma fluctuation structure
0122    info_dke_yp(2,['Vectorial form of the plasma fluctuations ',equil.id,'_',fluct.id,' is calculated.']);
0123    if C3POdisplay.fluctuations,testfitfluct_yp(equil_fit,fluct,fluct_fit);end
0124 end
0125 %
0126 % Calculation of Npar0 using Nphi0
0127 %
0128 % warning : this is valid only for m=0 cases
0129 %
0130 Bx_a0_fit = ppval(equil_fit.Bx_fit.pp_a0,rho0);
0131 Bx_an_fit = ppval(equil_fit.Bx_fit.pp_an,rho0);
0132 Bx_bn_fit = ppval(equil_fit.Bx_fit.pp_bn,rho0);
0133 %
0134 By_a0_fit = ppval(equil_fit.By_fit.pp_a0,rho0);
0135 By_an_fit = ppval(equil_fit.By_fit.pp_an,rho0);
0136 By_bn_fit = ppval(equil_fit.By_fit.pp_bn,rho0);
0137 %
0138 Bz_a0_fit = ppval(equil_fit.Bz_fit.pp_a0,rho0);
0139 Bz_an_fit = ppval(equil_fit.Bz_fit.pp_an,rho0);
0140 Bz_bn_fit = ppval(equil_fit.Bz_fit.pp_bn,rho0);
0141 %
0142 B_a0_fit = ppval(equil_fit.B_fit.pp_a0,rho0);
0143 B_an_fit = ppval(equil_fit.B_fit.pp_an,rho0);
0144 B_bn_fit = ppval(equil_fit.B_fit.pp_bn,rho0);
0145 %
0146 % Build interpolated magnetic fields
0147 %
0148 for jt = 1:length(theta0),
0149     [xBx] = calcval_yp(equil_fit,theta0(jt),Bx_a0_fit,Bx_an_fit,Bx_bn_fit);
0150     [xBy] = calcval_yp(equil_fit,theta0(jt),By_a0_fit,By_an_fit,By_bn_fit);
0151     [xBz] = calcval_yp(equil_fit,theta0(jt),Bz_a0_fit,Bz_an_fit,Bz_bn_fit);
0152     [xB] = calcval_yp(equil_fit,theta0(jt),B_a0_fit,B_an_fit,B_bn_fit);
0153     %
0154     xBzn = xBz./xB;
0155     %
0156     Npar0(jt) = Nphi0.*xBzn;
0157 end
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 toroidal position at launch
0163 rayinit.ym0 = m0;%Initial poloidal mode number
0164 rayinit.yn0 = n0;%Initial toroidal mode number
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 rayinit_theta = rayinit;
0178 rayinit_theta.ytheta0 = rayinit_theta.ytheta0(jt);
0179 rayinit_theta.yNpar0 = rayinit_theta.yNpar0(jt);
0180 %
0181 tstart = tic;
0182 wave_nofluct = main_C3PO_jd(dkepath,[id_wave,'_nofluct'],equil,equil_fit,rayinit_theta,waveparam,[],rayparam,C3POdisplay,C3POparam,[],[],0);clear mex;clear functions
0183 telapsed_ray_nofluct = toc(tstart);
0184 %
0185 info_dke_yp(2,'Ray trajectories calculated (interpolated magnetic equilibrium with no plasma fluctuations)');
0186 %
0187 tstart = tic;
0188 wave_fluct = main_C3PO_jd(dkepath,[id_wave,'_fluct'],equil,equil_fit,rayinit_theta,waveparam,[],rayparam,C3POdisplay,C3POparam,[],fluct_fit,0);
0189 telapsed_ray_fluct = toc(tstart);
0190 %
0191 info_dke_yp(2,'Ray trajectories calculated (interpolated magnetic equilibrium with plasma fluctuations)');
0192 %
0193 save_str = ['WAVE_',id_wave,'.mat'];
0194 save(save_str,'id_wave','wave_fluct','wave_nofluct','equil_fit','fluct');
0195 %
0196 info_dke_yp(2,'Wave parameters saved');
0197 %
0198 disp(['Elapsed time (s):',num2str(telapsed_ray_nofluct),' for ray calculations without density fluctuations (PCHIP interp. technique)']);
0199 disp(['Elapsed time (s):',num2str(telapsed_ray_fluct),' for ray calculations with density fluctuations (PCHIP interp. technique)']);
0200 %
0201 % --- display results ---
0202 %
0203 waves = {wave_nofluct,wave_fluct};
0204 rays = {wave_nofluct.rays{1},wave_fluct.rays{1}};
0205 %
0206 legs = {'No fluct.','With fluct.';...
0207        'nofluct','fluct'};
0208 %
0209 filename = ['Fig_',id_wave];
0210 %
0211 opt.p_opt = C3POdisplay.p_opt;
0212 opt.ntheta_fit = 65;
0213 opt.nrho_fit = 15;
0214 opt.propvar = 1;   
0215 %
0216 graph_comp_RT_jd(rays,legs,equil_fit,filename,opt)
0217 %
0218 diary4cvs_C3PO_yp(id_wave,dkepath,waves);% diary some results for CVS validation
0219 
0220 
0221 
0222 
0223 
0224 
0225 
0226 
0227 
0228 
0229 
0230 
0231