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 % 0019 % ***********************This part must be specified by the user, run make files if necessary) ***************************** 0020 % 0021 id_simul = 'LH_karney_test';%Simulation ID 0022 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0023 % 0024 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty) 0025 rho_S = [0.5];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty) 0026 % 0027 id_path = '';%For all paths used by DKE solver 0028 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0029 % 0030 id_equil = 'TScyl';%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 = 'UNIFORM10010020';%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 = 'NO_DISPLAY';%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 % 0057 %************************************************************************************************************************************ 0058 % 0059 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0 0060 %Option parameter for cross-comparison between old LH code: 0061 % - (1): 1/vpar dependence 0062 % - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd) 0063 wavestruct.opt_lh = 2; % [1,1] 0064 % 0065 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum, 0066 % or (Npar_lh,dNpar_lh) for Gaussian shape 0067 % 0068 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries 0069 % 0070 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh] 0071 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh] 0072 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh] 0073 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh] 0074 % 0075 % Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only 0076 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh] 0077 % 0078 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh] 0079 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0080 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0081 % 0082 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh] 0083 % (0) from local values Te and ne, (1) from central values Te0 and ne0 0084 % 0085 %************************************************************************************************************************************ 0086 % 0087 if exist('dmumpsmex');dkeparam.invproc = -2;disp('MUMPS mexfile is available -> matrix inversion uses this solver.');;end 0088 % 0089 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) 0090 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 0091 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution 0092 % 0093 %dkeparam.prec0_f = 1e-11;%to retrieve early LUKE 1.9.0 results 0094 %dkeparam.tn = 10000;%to retrieve early LUKE 1.9.0 results 0095 %dkeparam.dtn = 10000;%to retrieve early LUKE 1.9.0 results 0096 % 0097 dkeparam.psin_S = psin_S; 0098 dkeparam.rho_S = rho_S; 0099 % 0100 % RELATIVISTIC CASE 0101 % 0102 dkeparam.pnmax_S = 30; 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{1} = make_idealLHwave_jd(equil,wavestruct); 0108 % 0109 [Znorm_rel,Zcurr_rel,ZP0_rel] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0110 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource,dke_warnings] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);% for display with proc_luke_jd 0111 % 0112 % NON-RELATIVISTIC CASE 0113 % 0114 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants 0115 % 0116 betath = 0.001;%validated for NR limit 0117 equil.pTe = betath^2*mc2*ones(size(equil.pTe)); 0118 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi)); 0119 % 0120 wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0121 wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0122 % 0123 waves{1} = make_idealLHwave_jd(equil,wavestruct); 0124 % 0125 dkeparam.coll_mode = 0;%For comparison with theoretical results 0126 % 0127 dkeparam.nmhu_S = 101; 0128 dkeparam.np_S = 101; 0129 dkeparam.pnmax_S = 10; 0130 % 0131 [Znorm_nr_0_10,Zcurr_0_10,ZP0_0_10] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0132 % 0133 dkeparam.nmhu_S = 201; 0134 dkeparam.np_S = 201; 0135 dkeparam.pnmax_S = 20; 0136 % 0137 [Znorm_nr_0_20,Zcurr_0_20,ZP0_0_20] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0138 % 0139 dkeparam.coll_mode = 2;%For comparison with theoretical results 0140 % 0141 [Znorm_nr_2_20,Zcurr_2_20,ZP0_2_20] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0142 % 0143 j_r_k = '0.003732'; 0144 p_r_k = '0.0001256'; 0145 eta_r_k = '29.72'; 0146 % 0147 j_nr_0_10_k = '0.05754'; 0148 p_nr_0_10_k = '0.004011'; 0149 eta_nr_0_10_k = '14.34'; 0150 % 0151 j_nr_0_20_k = '0.05759'; 0152 p_nr_0_20_k = '0.004012'; 0153 eta_nr_0_20_k = '14.35'; 0154 % 0155 j_nr_2_20_k = '0.07092'; 0156 p_nr_2_20_k = '0.004294'; 0157 eta_nr_2_20_k = '16.52'; 0158 % 0159 %************************************************************************************************************************************ 0160 % 0161 format 0162 % 0163 delete res_karney 0164 % 0165 diary res_karney 0166 % 0167 disp(['Comparison LUKE/Karney with tn = ',int2str(dkeparam.tn),', dtn = ',int2str(dkeparam.dtn),' and nit_f = ',int2str(dkeparam.nit_f)]) 0168 disp(['----------------------']) 0169 % 0170 disp(['Relativistic case (4,7,1) : j = ',num2str(Zcurr_rel.x_0),'/',j_r_k,' ; P = ',num2str(ZP0_rel.x_rf_fsav),'/',p_r_k,' ; j/P = ',num2str(Zcurr_rel.x_0./ZP0_rel.x_rf_fsav),'/',eta_r_k]) 0171 disp(['Non Relativistic case (3,5,1) pmax = 10, coll_mode = 0 : j = ',num2str(Zcurr_0_10.x_0),'/',j_nr_0_10_k,' ; P = ',num2str(ZP0_0_10.x_rf_fsav),'/',p_nr_0_10_k,' ; j/P = ',num2str(Zcurr_0_10.x_0./ZP0_0_10.x_rf_fsav),'/',eta_nr_0_10_k]) 0172 disp(['Non Relativistic case (3,5,1) pmax = 20, coll_mode = 0 : j = ',num2str(Zcurr_0_20.x_0),'/',j_nr_0_20_k,' ; P = ',num2str(ZP0_0_20.x_rf_fsav),'/',p_nr_0_20_k,' ; j/P = ',num2str(Zcurr_0_20.x_0./ZP0_0_20.x_rf_fsav),'/',eta_nr_0_20_k]) 0173 disp(['Non Relativistic case (3,5,1) pmax = 20, coll_mode = 2 : j = ',num2str(Zcurr_2_20.x_0),'/',j_nr_2_20_k,' ; P = ',num2str(ZP0_2_20.x_rf_fsav),'/',p_nr_2_20_k,' ; j/P = ',num2str(Zcurr_2_20.x_0./ZP0_2_20.x_rf_fsav),'/',eta_nr_2_20_k]) 0174 % 0175 diary off 0176 % 0177 %************************************************************************************************************************************ 0178 % 0179 str = [path_simul,'LUKE_RESULTS_',id_equil,'_',id_simul,'.mat']; 0180 eval(['save ',str]); 0181 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']); 0182 % 0183 proc_luke_jd(load(str),1,2,0.1,20);% for the relativist case only