rundke_val

Script for running the DKE solver (can be modified by the user for specific simulations)
by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)

0001 %Script for running the DKE solver (can be modified by the user for specific simulations)
0002 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)
0003 %
0004 clear all
0005 clear mex
0006 clear functions
0007 close all
0008 warning off
0009 global nfig
0010 %
0011 permission = test_permissions_yp;
0012 %
0013 if ~permission 
0014     disp('Please move the script to a local folder where you have write permission before to run it')
0015     return;
0016 end
0017 %
0018 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0019 %
0020 id_simul = 'LH_bounce_val';%Simulation ID
0021 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0022 %
0023 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0024 rho_S = [];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0025 %
0026 id_path = '';%For all paths used by DKE solver
0027 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0028 %
0029 id_equil_cyl = 'TScyl';%For plasma equilibrium
0030 id_equil = 'TScirc_e1';%For plasma equilibrium
0031 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0032 %
0033 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0034 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0035 %
0036 id_display = 'PARTIAL_VISUAL';%For output code display
0037 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0038 %
0039 id_ohm = '';%For Ohmic electric contribution
0040 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0041 %
0042 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0043 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0044 %
0045 id_transpfaste = '';%For fast electron radial transport
0046 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0047 %
0048 id_ripple = '';%For fast electron magnetic ripple losses
0049 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0050 %
0051 %************************************************************************************************************************************
0052 %************************************************************************************************************************************
0053 %************************************************************************************************************************************
0054 %
0055 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple);
0056 [equil_cyl] = load_structures_yp('equil',id_equil_cyl,path_equil);
0057 %
0058 %************************************************************************************************************************************
0059 %
0060 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0061 %Option parameter for cross-comparison between old LH code:
0062 %    - (1): 1/vpar dependence
0063 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0064 wavestruct.opt_lh = 2; % [1,1]
0065 %
0066 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0067 % or (Npar_lh,dNpar_lh) for Gaussian shape
0068 %
0069 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0070 %
0071 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0072 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0073 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0074 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0075 %
0076 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0077 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0078 %
0079 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0080 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0081 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0082 %
0083 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0084 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0085 %
0086 %************************************************************************************************************************************
0087 %
0088 if exist('dmumpsmex');dkeparam.invproc = -2;end
0089 %
0090 dkeparam.boundary_mode_f = 0;%Number of points where the Maxwellian distribution is enforced from p = 0 (p=0, free conservative mode but param_inv(1) must be less than 1e-4, otherwise 1e-3 is OK most of the time. Sensitive to the number of points in p)
0091 dkeparam.norm_mode_f = 1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0092 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0093 %
0094 dkeparam.psin_S = psin_S;
0095 %
0096 iar_list = [0.1,0.2,0.3,0.5,0.8];
0097 %
0098 % RELATIVISTIC CASE
0099 %
0100 dkeparam.pnmax_S = 30;
0101 dkeparam.np_S = 201;
0102 dkeparam.nmhu_S = 201;
0103 %
0104 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0105 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0106 %
0107 waves_cyl{1} = make_idealLHwave_jd(equil_cyl,wavestruct);
0108 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0109 %
0110 ap = equil.ptx(end,1);
0111 rho_list = iar_list*equil.Rp/ap;
0112 %
0113 nrho = length(rho_list);
0114 j_rel = NaN(1,nrho);
0115 p_rel = NaN(1,nrho);
0116 j_rel_cyl = NaN(1,nrho);
0117 p_rel_cyl = NaN(1,nrho);
0118 %
0119 for irho = 1:nrho,
0120     %
0121     dkeparam.bounce_mode = 1;
0122     dkeparam.mhu_S = [];
0123     %
0124     dkeparam.rho_S = rho_list(irho);
0125     %
0126     [Znorm_rel,Zcurr_rel,ZP0_rel,dke_out_rel,radialDKE_rel,equilDKE_rel,momentumDKE_rel] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0127     %
0128     if dke_out_rel.residu_f{end}(end) <= dkeparam.prec0_f,
0129         j_rel(irho) = Zcurr_rel.x_0_vfsav;
0130         p_rel(irho) = ZP0_rel.x_rf_fsav.';
0131     end  
0132     %
0133     dkeparam.bounce_mode = 0;
0134     dkeparam.mhu_S = momentumDKE_rel.mhu_S;
0135     %
0136     dkeparam.rho_S = 0.5;
0137     %
0138     [Znorm_rel_cyl,Zcurr_rel_cyl,ZP0_rel_cyl,dke_out_rel_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm,waves_cyl,transpfaste,ripple,[],[]);
0139     %
0140     if dke_out_rel_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0141         j_rel_cyl(irho) = Zcurr_rel_cyl.x_0;
0142         p_rel_cyl(irho) = ZP0_rel_cyl.x_rf_fsav.';
0143     end  
0144     %
0145 end
0146 %
0147 % NON-RELATIVISTIC CASE
0148 %
0149 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0150 %
0151 betath = 0.001;%validated for NR limit
0152 %
0153 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0154 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0155 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0156 %
0157 equil_cyl.pTe = betath^2*mc2*ones(size(equil_cyl.pTe));
0158 equil_cyl.pzTi = betath^2*mc2*ones(size(equil_cyl.pzTi));
0159 %equil_cyl.pzTi = 1e-10*ones(size(equil_cyl.pzTi));
0160 %
0161 wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0162 wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0163 %
0164 waves_cyl{1} = make_idealLHwave_jd(equil_cyl,wavestruct);
0165 waves{1} = make_idealLHwave_jd(equil,wavestruct);
0166 %
0167 dkeparam.pnmax_S = 10;
0168 %
0169 j_nr = NaN(1,nrho);
0170 p_nr = NaN(1,nrho);
0171 j_nr_cyl = NaN(1,nrho);
0172 p_nr_cyl = NaN(1,nrho);
0173 %
0174 for irho = 1:nrho,
0175     %
0176     dkeparam.bounce_mode = 1;
0177     dkeparam.mhu_S = [];
0178     %
0179     dkeparam.rho_S = rho_list(irho);
0180     %
0181     [Znorm_nr,Zcurr_nr,ZP0_nr,dke_out_nr,radialDKE_nr,equilDKE_nr,momentumDKE_nr] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0182     %
0183     if dke_out_nr.residu_f{end}(end) <= dkeparam.prec0_f,
0184         j_nr(irho) = Zcurr_nr.x_0_vfsav;
0185         p_nr(irho) = ZP0_nr.x_rf_fsav.';
0186     end  
0187     %
0188     dkeparam.bounce_mode = 0;
0189     dkeparam.mhu_S = momentumDKE_nr.mhu_S;
0190     %
0191     dkeparam.rho_S = 0.5;
0192     %
0193     [Znorm_nr_cyl,Zcurr_nr_cyl,ZP0_nr_cyl,dke_out_nr_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm,waves_cyl,transpfaste,ripple,[],[]);
0194     %
0195     if dke_out_nr_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0196         j_nr_cyl(irho) = Zcurr_nr_cyl.x_0;
0197         p_nr_cyl(irho) = ZP0_nr_cyl.x_rf_fsav.';
0198     end  
0199     %
0200 end
0201 %
0202 eta_rel_cyl = j_rel_cyl./p_rel_cyl;
0203 eta_rel = j_rel./p_rel;
0204 %
0205 eta_nr_cyl = j_nr_cyl./p_nr_cyl;
0206 eta_nr = j_nr./p_nr;
0207 %
0208 %************************************************************************************************************************************
0209 %
0210 format
0211 %
0212 delete res_LH_bounce_val
0213 %
0214 diary res_LH_bounce_val
0215 %
0216 disp(['Case 1 (NR)'])
0217 disp(['----------------------'])
0218 disp(['--> eps = 0.0 : eta/eta_cyl = ',num2str(eta_nr_cyl)]) 
0219 disp(['--> eps = 0.1 : eta/eta_cyl = ',num2str(eta_nr(1))]) 
0220 disp(['--> eps = 0.2 : eta/eta_cyl = ',num2str(eta_nr(2))]) 
0221 disp(['--> eps = 0.3 : eta/eta_cyl = ',num2str(eta_nr(3))]) 
0222 disp(['--> eps = 0.5 : eta/eta_cyl = ',num2str(eta_nr(4))]) 
0223 disp(['--> eps = 0.8 : eta/eta_cyl = ',num2str(eta_nr(5))]) 
0224 %
0225 %
0226 disp(' ')
0227 disp(['Case 2 (FR)'])
0228 disp(['----------------------'])
0229 disp(['--> eps = 0.0 : eta/eta_cyl = ',num2str(eta_rel_cyl)]) 
0230 disp(['--> eps = 0.1 : eta/eta_cyl = ',num2str(eta_rel(1))]) 
0231 disp(['--> eps = 0.2 : eta/eta_cyl = ',num2str(eta_rel(2))]) 
0232 disp(['--> eps = 0.3 : eta/eta_cyl = ',num2str(eta_rel(3))]) 
0233 disp(['--> eps = 0.5 : eta/eta_cyl = ',num2str(eta_rel(4))]) 
0234 disp(['--> eps = 0.8 : eta/eta_cyl = ',num2str(eta_rel(5))]) 
0235 %
0236 diary off
0237 %
0238 %************************************************************************************************************************************
0239 %
0240 %
0241 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0242 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);