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 opt_proc = 0;%test for numerical performance test with number of processors 0012 Dr0 = 1; 0013 % 0014 p_opt = 2; 0015 % 0016 permission = test_permissions_yp; 0017 % 0018 if ~permission 0019 disp('Please move the script to a local folder where you have write permission before to run it') 0020 return; 0021 end 0022 % 0023 % ***********************This part must be specified by the user, run make files if necessary) ***************************** 0024 % 0025 id_simul = ['LH_karney_3D_splitting_',num2str(Dr0)];%Simulation ID 0026 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0027 % 0028 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty) 0029 rho_S = [0:0.05:0.65,0.7:0.1:1];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty) 0030 % 0031 id_path = '';%For all paths used by DKE solver 0032 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0033 % 0034 id_equil = 'TScyl';%For plasma equilibrium 0035 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0036 % 0037 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters 0038 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0039 % 0040 id_display = 'PARTIAL_VISUAL';%For output code display 0041 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0042 % 0043 id_ohm = '';%For Ohmic electric contribution 0044 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0045 % 0046 ids_wave = {''};%For RF waves contribution (put all the type of waves needed) 0047 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0048 % 0049 id_transpfaste = 'nomagneticturbulence';%For fast electron radial transport 0050 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0051 % 0052 id_ripple = '';%For fast electron magnetic ripple losses 0053 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0054 % 0055 %************************************************************************************************************************************ 0056 %************************************************************************************************************************************ 0057 %************************************************************************************************************************************ 0058 % 0059 [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); 0060 % 0061 %************************************************************************************************************************************ 0062 % 0063 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0 0064 %Option parameter for cross-comparison between old LH code: 0065 % - (1): 1/vpar dependence 0066 % - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd) 0067 wavestruct.opt_lh = 2; % [1,1] 0068 % 0069 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum, 0070 % or (Npar_lh,dNpar_lh) for Gaussian shape 0071 % 0072 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries 0073 % 0074 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0075 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0076 % 0077 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh] 0078 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh] 0079 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh] 0080 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh] 0081 % 0082 % Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only 0083 wavestruct.yD0_in_c_lh = [2];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh] 0084 % 0085 wavestruct.yD0_in_lh_prof = [1];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh] 0086 wavestruct.ypeak_lh = [0.4];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0087 wavestruct.ywidth_lh = [0.1];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0088 % 0089 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh] 0090 % (0) from local values Te and ne, (1) from central values Te0 and ne0 0091 % 0092 %************************************************************************************************************************************ 0093 % 0094 if exist('dmumpsmex');dkeparam.invproc = -2;end 0095 % 0096 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) 0097 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 0098 dkeparam.timevol = 1;%to calculate moments at all internal times 0099 % 0100 dkeparam.psin_S = psin_S; 0101 dkeparam.rho_S = rho_S; 0102 % 0103 %betath = 0.001; 0104 % 0105 %[qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants 0106 % 0107 %equil.pTe = betath^2*mc2*ones(size(equil.pTe)); 0108 %equil.pzTi = betath^2*mc2*ones(size(equil.pzTi)); 0109 %equil.pzTi = 1e-10*ones(size(equil.pzTi)); 0110 % 0111 waves{1} = make_idealLHwave_jd(equil,wavestruct); 0112 % 0113 %dkedisplay.display_mode = 0;%If -1, only the radius at which the residu is the largest is displayed 0114 % 0115 tfinal = 2000; 0116 % 0117 dtn_list = [1,10,100,1000]; 0118 dtn_list = [10,100,1000]; 0119 coll_mode_list = [0,2]; 0120 optsplit_list = [0,1]; 0121 % 0122 ncoll = length(coll_mode_list); 0123 nopt = length(optsplit_list); 0124 ndtn = length(dtn_list); 0125 % 0126 dkeparam.tn = tfinal; 0127 transpfaste.Dr0 = Dr0; 0128 % 0129 I = cell(ndtn,nopt,ncoll); 0130 % 0131 for icoll = 1:ncoll, 0132 for iopt = 1:nopt, 0133 for idtn = ndtn:-1:1, 0134 % 0135 dkeparam.coll_mode = coll_mode_list(icoll); 0136 dkeparam.opsplit = optsplit_list(iopt); 0137 dkeparam.dtn = dtn_list(idtn); 0138 % 0139 [dummy,Zcurr,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0140 % 0141 nit = tfinal/dkeparam.dtn; 0142 I{idtn,iopt,icoll} = NaN(1,nit); 0143 % 0144 for it = 1:nit, 0145 if dkeparam.coll_mode ~= 2 || dke_out.residu_f{it}(end) <= dkeparam.prec0_f, 0146 I{idtn,iopt,icoll}(it) = Zcurr(it).I_tot; 0147 end 0148 end 0149 % 0150 disp(['- collmode: ',int2str(dkeparam.coll_mode),',- opsplit: ',int2str(dkeparam.opsplit),',dtn: ',int2str(dtn_list(idtn)),' done.']) 0151 end 0152 end 0153 end 0154 % 0155 %************************************************************************************************************************************ 0156 % 0157 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']); 0158 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']); 0159 % 0160 proc_rundke_splitting