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_simul0 = 'LH_karney';%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 = [0.05:0.1:0.9];%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 = 'TScirc_peak';%For plasma equilibrium 0030 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0031 % 0032 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters 0033 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0034 % 0035 id_display = 'PARTIAL_VISUAL';%For output code display 0036 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0037 % 0038 id_ohm = '';%For Ohmic electric contribution 0039 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0040 % 0041 ids_wave = {''};%For RF waves contribution (put all the type of waves needed) 0042 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0043 % 0044 id_transpfaste = '';%For fast electron radial transport 0045 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0046 % 0047 id_ripple = 'TSspecific';%For fast electron magnetic ripple losses 0048 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path 0049 % 0050 %************************************************************************************************************************************ 0051 %************************************************************************************************************************************ 0052 %************************************************************************************************************************************ 0053 % 0054 [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); 0055 % 0056 %************************************************************************************************************************************ 0057 % 0058 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0 0059 %Option parameter for cross-comparison between old LH code: 0060 % - (1): 1/vpar dependence 0061 % - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd) 0062 wavestruct.opt_lh = 2; % [1,1] 0063 % 0064 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum, 0065 % or (Npar_lh,dNpar_lh) for Gaussian shape 0066 % 0067 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries 0068 % 0069 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0070 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh] 0071 % 0072 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh] 0073 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh] 0074 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh] 0075 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh] 0076 % 0077 % Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only 0078 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh] 0079 % 0080 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh] 0081 wavestruct.ypeak_lh = [0.4];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0082 wavestruct.ywidth_lh = [0.2];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh] 0083 % 0084 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh] 0085 % (0) from local values Te and ne, (1) from central values Te0 and ne0 0086 % 0087 %************************************************************************************************************************************ 0088 % 0089 if exist('dmumpsmex');dkeparam.invproc = -2;end 0090 % 0091 dkeparam.boundary_mode_f = 1;%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) 0092 dkeparam.norm_mode_f = 0;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes 0093 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution 0094 % 0095 dkeparam.np_S = 201; 0096 dkeparam.nmhu_S = 201; 0097 dkeparam.pnmax_S = 30; 0098 % 0099 dkeparam.psin_S = psin_S; 0100 dkeparam.rho_S = rho_S; 0101 % 0102 waves{1} = make_idealLHwave_jd(equil,wavestruct); 0103 % 0104 % case with no ripple 0105 % 0106 id_simul_off = [id_simul0,'_rippleoff']; 0107 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul_off,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,[],[],[]); 0108 str = ['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul_off,'.mat Znorm Zcurr ZP0 dke_out radialDKE equilDKE momentumDKE gridDKE Zmomcoef Zbouncecoef Zmripple mksa']; 0109 eval(str); 0110 % 0111 % case with ripple 0112 % 0113 id_simul_on = [id_simul0,'_rippleon']; 0114 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul_on,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]); 0115 str = ['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul_on,'.mat Znorm Zcurr ZP0 dke_out radialDKE equilDKE momentumDKE gridDKE Zmomcoef Zbouncecoef Zmripple mksa']; 0116 eval(str); 0117 % 0118 %************************************************************************************************************************************ 0119 % 0120 vparmin = 1/(waves{1}.rays{1}.sNpar(1)+imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref; 0121 vparmax = 1/(waves{1}.rays{1}.sNpar(1)-imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref; 0122 subtitle = ['Tore Supra, v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{||LH} = ',num2str(wavestruct.yD0_in_c_lh),', Zi = ',num2str(max(equilDKE.xZeff))]; 0123 % 0124 % 0125 print_opt = 2; 0126 proc_luke_jd(load(['DKE_RESULTS_',id_equil,'_',id_simul_off]),find(dke_out.xMRR_flux == max(dke_out.xMRR_flux)),print_opt,0.1,20,subtitle);%Display where magnetic ripple losses peak 0127 proc_luke_jd(load(['DKE_RESULTS_',id_equil,'_',id_simul_on]),find(dke_out.xMRR_flux == max(dke_out.xMRR_flux)),print_opt,0.1,20,subtitle);%Display where magnetic ripple losses peak 0128 %