make_wave_JETlikeflucttest_pureripple

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