rundke_cluster

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_karney_cluster';%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.5];%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 = 'TScyl';%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 = 'UNIFORM10010020';%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 = 'NO_DISPLAY';%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 = '';%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.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0070 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0071 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0072 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0073 %
0074 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0075 wavestruct.yD0_in_c_lh = [1];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0076 %
0077 wavestruct.yD0_in_lh_prof = [0];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0078 wavestruct.ypeak_lh = [NaN];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0079 wavestruct.ywidth_lh = [NaN];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0080 %
0081 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0082 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0083 %
0084 %************************************************************************************************************************************
0085 %
0086 if exist('dmumpsmex');dkeparam.invproc = -2;end
0087 %
0088 % Distributed computing parameters (scheduler choice)
0089 %
0090 dkeparam.clustermode.coll_dke_jd.scheduler = 'local';%'' or 'local', 'pbspro' or otherwise
0091 %
0092 dkeparam.clustermode.eecoll_dke_yp.scheduler = 'local';%'' or 'local', 'pbspro' or otherwise
0093 %
0094 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)
0095 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
0096 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0097 %
0098 dkeparam.psin_S = psin_S;
0099 dkeparam.rho_S = rho_S;
0100 %
0101 j_r_k = '0.003732';
0102 p_r_k = '0.0001256';
0103 eta_r_k = '29.72';
0104 %
0105 j_nr_0_10_k = '0.05754';
0106 p_nr_0_10_k = '0.004011';
0107 eta_nr_0_10_k = '14.34';
0108 %
0109 j_nr_0_20_k = '0.05759';
0110 p_nr_0_20_k = '0.004012';
0111 eta_nr_0_20_k = '14.35';
0112 %
0113 j_nr_2_20_k = '0.07092';
0114 p_nr_2_20_k = '0.004294';
0115 eta_nr_2_20_k = '16.52';
0116 %
0117 iclustermode_min = 1;%1: sequential, 2: MDCE, 3: LDCE, 4: PARFOR, 5: GPU
0118 iclustermode_max = 4;%1: sequential, 2: MDCE, 3: LDCE, 4: PARFOR, 5: GPU
0119 %
0120 for iclustermode = iclustermode_min:iclustermode_max,
0121     %
0122     dkeparam.clustermode.coll_dke_jd.scheduler.mode = iclustermode-1;%0: sequential, 1: MDCE, 2: LDCE, 3: PARFOR, 4: GPU
0123     dkeparam.clustermode.coll_dke_jd.scheduler.memory = 500;%memory allocated to batch runs, in MB
0124     %
0125     dkeparam.clustermode.eecoll_dke_yp.scheduler.mode = iclustermode-1;%0: sequential, 1: MDCE, 2: LDCE, 3: PARFOR, 4: GPU
0126     dkeparam.clustermode.eecoll_dke_yp.scheduler.memory = 500;%memory allocated to batch runs, in MB
0127     %
0128     % RELATIVISTIC CASE
0129     %
0130     dkeparam.pnmax_S = 30;
0131     %
0132     wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0133     wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0134     %
0135     waves{1} = make_idealLHwave_jd(equil,wavestruct);
0136     %
0137     [Znorm_rel{iclustermode},Zcurr_rel{iclustermode},ZP0_rel{iclustermode}] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0138     %
0139     % NON-RELATIVISTIC CASE
0140     %
0141     [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0142     %
0143     betath = 0.001;%validated for NR limit
0144     equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0145     equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0146     %
0147     wavestruct.yvparmin_lh = [3];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0148     wavestruct.yvparmax_lh = [5];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0149     %
0150     waves{1} = make_idealLHwave_jd(equil,wavestruct);
0151     %
0152     dkeparam.coll_mode = 0;%For comparison with theoretical results
0153     %
0154     dkeparam.nmhu_S = 101;
0155     dkeparam.np_S = 101;
0156     dkeparam.pnmax_S = 10;
0157     %
0158     [Znorm_nr_0_10{iclustermode},Zcurr_0_10{iclustermode},ZP0_0_10{iclustermode}] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0159     %
0160     dkeparam.nmhu_S = 201;
0161     dkeparam.np_S = 201;
0162     dkeparam.pnmax_S = 20;
0163     %
0164     [Znorm_nr_0_20{iclustermode},Zcurr_0_20{iclustermode},ZP0_0_20{iclustermode}] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0165     %
0166     dkeparam.coll_mode = 2;%For comparison with theoretical results
0167     %
0168     [Znorm_nr_2_20{iclustermode},Zcurr_2_20{iclustermode},ZP0_2_20{iclustermode}] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0169 end
0170 %
0171 clc
0172 %
0173 disp(['--------------------------------------------------------------------------------------------'])
0174 %
0175 format
0176 %
0177 delete res_karney_cluster
0178 %
0179 diary res_karney_cluster
0180 %
0181 for iclustermode = iclustermode_min:iclustermode_max,  
0182     disp(' ')
0183     if iclustermode == 1,
0184         disp(['Comparison LUKE/Karney, sequential calculations'])
0185     elseif iclustermode == 2,
0186         disp(['Comparison LUKE/Karney, Matlab Distributed Computing Toolbox calculations (MDCE)'])
0187     elseif iclustermode == 3,
0188         disp(['Comparison LUKE/Karney, LUKE Distributed Computing Toolbox calculations (LDCE)'])
0189     elseif iclustermode == 4,
0190         disp(['Comparison LUKE/Karney, with PARFOR loop calculations'])
0191     elseif iclustermode == 5,
0192         disp(['Comparison LUKE/Karney, with GPU calculations'])
0193     end
0194     %
0195     disp(['----------------------------------------------'])
0196     %
0197     disp(['Relativistic case (4,7,1) : j = ',num2str(Zcurr_rel{iclustermode}.x_0),'/',j_r_k,' ; P = ',num2str(ZP0_rel{iclustermode}.x_rf_fsav),'/',p_r_k,' ; j/P = ',num2str(Zcurr_rel{iclustermode}.x_0./ZP0_rel{iclustermode}.x_rf_fsav),'/',eta_r_k]) 
0198     disp(['Non Relativistic case (3,5,1) pmax = 10, coll_mode = 0 : j = ',num2str(Zcurr_0_10{iclustermode}.x_0),'/',j_nr_0_10_k,' ; P = ',num2str(ZP0_0_10{iclustermode}.x_rf_fsav),'/',p_nr_0_10_k,' ; j/P = ',num2str(Zcurr_0_10{iclustermode}.x_0./ZP0_0_10{iclustermode}.x_rf_fsav),'/',eta_nr_0_10_k]) 
0199     disp(['Non Relativistic case (3,5,1) pmax = 20, coll_mode = 0 : j = ',num2str(Zcurr_0_20{iclustermode}.x_0),'/',j_nr_0_20_k,' ; P = ',num2str(ZP0_0_20{iclustermode}.x_rf_fsav),'/',p_nr_0_20_k,' ; j/P = ',num2str(Zcurr_0_20{iclustermode}.x_0./ZP0_0_20{iclustermode}.x_rf_fsav),'/',eta_nr_0_20_k]) 
0200     disp(['Non Relativistic case (3,5,1) pmax = 20, coll_mode = 2 : j = ',num2str(Zcurr_2_20{iclustermode}.x_0),'/',j_nr_2_20_k,' ; P = ',num2str(ZP0_2_20{iclustermode}.x_rf_fsav),'/',p_nr_2_20_k,' ; j/P = ',num2str(Zcurr_2_20{iclustermode}.x_0./ZP0_2_20{iclustermode}.x_rf_fsav),'/',eta_nr_2_20_k]) 
0201 end
0202 %
0203 diary off
0204 %
0205 %************************************************************************************************************************************
0206 %
0207 disp(' ')
0208 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0209 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);