Create the strucure dkeparam for CRONOS by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu) --------- Load first default parameters ------------
0001 function [dkeparam] = make_dkeparam_cronos(par) 0002 % 0003 % Create the strucure dkeparam for CRONOS 0004 % 0005 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu) 0006 % 0007 % --------- Load first default parameters ------------ 0008 % 0009 if strcmp(par.luke_mode,'lh') || strcmp(par.luke_mode,'mixte');%For the Lower Hybrid wave simulations 0010 % 0011 %Mesh size calculation 0012 % 0013 dkeparam.id = 'cronos_LH'; 0014 % 0015 dkeparam.grid_uniformity = 0;%Pitch-angle and momentum grids uniformity: (0) no, (1) yes (-1 for specific test for DKE, transient and pitch-angle grid only) 0016 dkeparam.grid_method = 1;%Pitch-angle grid calculation technique: (0): old method, (1) new method 0017 dkeparam.nmhu_S = 121;%Number of angular values for the flux grid (WARNING: odd number) (>100 recommended for accurate run) 0018 dkeparam.np_S = 151;%Number of momentum values for the flux grid (pnmax/np < 0.1 highly recommended for accurate run) WARNING: must be adjusted so that effective momentum point number is larger than thepitch-angle one. 0019 dkeparam.pnmax_S = 30;%Maximum value of the flux momentum grid (pth_ref) 0020 dkeparam.rho_S = 20;% number of points for radial grid calculation 0021 dkeparam.rho0 = 0.01;%option for extra grid point near 0 (value above 0, none if 0) 0022 % 0023 if isfield(par,'eco') && strcmp(par.eco,'Yes') 0024 disp('memory economic mode on') 0025 dkeparam.nmhu_S = 101;%Number of angular values for the flux grid (WARNING: odd number) (>100 recommended for accurate run) 0026 dkeparam.np_S = 121;%Number of momentum values for the flux grid (pnmax/np < 0.1 highly recommended for accurate run) WARNING: must be adjusted so that effective momentum point number is larger than thepitch-angle one. 0027 dkeparam.pnmax_S = 20;%Maximum value of the flux momentum grid (pth_ref) 0028 par.rho_S = '15'; 0029 end 0030 % 0031 %Pitch-angle grid 0032 % 0033 dkeparam.sfac0 = 1; % reduction factor for fine grid near mhu = 0 0034 dkeparam.sfac1 = 5; % reduction factor for fine grid near abs(mhu) = 1 0035 dkeparam.sfact = 1; % reduction factor for fine grid near abs(mhu) = mhut 0036 % 0037 dkeparam.snmhu0 = 5; % number of fine grid points near mhu = 0 0038 dkeparam.snmhu1 = 1; % number of fine grid points near abs(mhu) = 1 0039 dkeparam.snmhut = 1; % number of fine grid points near abs(mhu) = mhut 0040 % 0041 dkeparam.grid_power = 1.5; % power for mhu grid step varations (grid_power > 1) (Note: for grid_method = 1 only) 0042 % 0043 %Momentum grid 0044 % 0045 dkeparam.np_S_ref = 20;%position in the mpomentum grid where dpn_S is twice the minimum value (must be ajusted with snp parameter) 0046 dkeparam.snp = 1;%pgridfactor (Uniform grid: pgridfactor = 1, non-uniform grid: pgridfactor > 1) 0047 dkeparam.spfac = 20;%reduction factor for the momentum supergrid used in collision integrals 0048 % 0049 %Normalized radii at which the distribution is calculated (DKE or FP) 0050 %The syntax rho='all' may be also used if all mhu values are considered. 0051 % 0052 dkeparam.dpsin_dke = 0.01;%maximum difference between two consecutive positions used for radial derivatives 0053 % 0054 dkeparam.ref_mode = 0;%Reference values: (0) from core Te and ne, (1) from Te and ne at the normalized radius (paramater only active for a single radial value, otherwise ref_mode = 0) 0055 % 0056 % Dql and collision operator 0057 % 0058 dkeparam.coll_mode = 2;%Relativistic collision operator: (0) Relativistic Maxwellian background, (1) High-velocity limit, (2) Linearized Belaiev-Budker (3) Non-relativistic Lorentz model, (4) Non-relativistic Maxwellian background 0059 % 0060 % Ray-tracing and RF waves 0061 % 0062 dkeparam.rt_mode = 4;%RT-DKE power consistency method 0063 dkeparam.abs_lim = 0.95;%minimum required absorption fraction 0064 dkeparam.nk_abs = 5;%maximum number of iterations on a_sdNpar for full absorption 0065 dkeparam.nit_rf = 40;%Maximum number of iteration in ray-tracing. Works only if rt_mode = 1 and more than one radial point (default = 1) 0066 dkeparam.prec_rf = 1e-2;%relative precision in RT-DKE iterations 0067 %dkeparam.n_rf_list = 0;%List of harmonics in the Dql calculations 0068 dkeparam.delta_rt = 0.5;%smoothering parameter of quasilinear iterations 0069 dkeparam.Dparparmax = 20;%Integration time step (1/nhuth_ref) (dtn > 0: fully implicit time differencing scheme , dtn < 0: Crank-Nicholson time differencing scheme) (WARNING: dtn less than 1 when fast electron radial transport is considered) 0070 dkeparam.npresmin = 1;%minimum number of resonant point per ray 0071 % 0072 % Run mode 0073 % 0074 dkeparam.bounce_mode = 1;%Bounce averaging option: (0) no bounce averaging, (1) bounce averaging full implicit method with 15 diagonals 0075 dkeparam.dke_mode = 0;%Fokker-Planck equation (0), Drift Kinetic equation (1) 0076 dkeparam.equil_mode = 1;%Equilibrium option: (0) circular concentric flux-surfaces, (1) numerical equilibrium 0077 dkeparam.syn_mode = 0;%Synergy calculation: (0) none, (1) yes: the code runs twice, including once without RF or electric fields. (only if dke_mode = 1) 0078 % 0079 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 0080 dkeparam.boundary_mode_f = 2;%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) 0081 dkeparam.prec0_f = 1e-11;%Convergence level of the residu 0082 dkeparam.nit_f = 50;%Maximum number of iterations for f (> 1) 0083 dkeparam.initguess_f = 0;%Initial guess for the distribution function: (0) Maxwellian, (1) Fuchs model with Tperp for LH current drive only, (2) from a file 0084 dkeparam.initfile_f = '';%If initguess_f = 2. The file must be in the directory Project_DKE, and created by dke_4_1yp.m 0085 % 0086 dkeparam.boundary_mode_g = 0;%Number of points where the first distribution function for the Lorentz model 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) 0087 dkeparam.prec0_g = 1e-10;%Convergence level of the residu (if NaN, iterations stop with nit_g) 0088 dkeparam.nit_g = 100;%Maximum number of iterations for g 0089 % 0090 dkeparam.dtn = NaN;%Integration time step (1/nhuth_ref) (dtn > 0: fully implicit time differencing scheme , dtn < 0: Crank-Nicholson time differencing scheme) (WARNING: dtn less than 1 when fast electron radial transport is considered) 0091 dkeparam.opsplit = 0;%Operator splitting between momentum and radial dynamics (for tests only) (opsplit = 0: no operator splitting , opsplit = 1: operator splitting). If this field does not exist, it is equivalent to dkeparam.opsplit = 0. 0092 % 0093 % (3) Matlab built-in qmr iterative method (MatLab built-in partial LU factorization used) 0094 % (2) Matlab built-in bicg iterative method (MatLab built-in partial LU factorization used) 0095 % (1) Matlab built-in cgs iterative method (MatLab built-in partial LU factorizsation used) 0096 % (0) Matlab built-in direct inversion without any LU factorization (NOT RECOMMANDED !) 0097 % (-1): MUMPS parallel (or sequential) fortran 95 direct inversion (LU matrix factorization done for each time a matrix inversion is called) 0098 % (-2): MUMPSMEX sequential fortran 95 inversion using MUMPS LU matrix factorization 0099 % (-3): SUPERLUMEX sequential direct matrix inversion (LU matrix factorization done for each time a matrix inversion is called) 0100 % (-30) PETSc interface (see data structure below) 0101 dkeparam.invproc = -2;%Iterative matrix inversion methods 0102 dkeparam.ludroptol = 1e-5;%Drop tolerance for the built-in MatLab partial LU decomposition (recipes: droptol = max(1/dtn/nr,param_inv(1))) 0103 dkeparam.mlmaxitinv = 4;%Maximum number of iterations for iterative built-in MatLab matrix inversion method (between 1 and 20 usualy, if empty, mlmaxitinv = min(matrixdim,20). See matlab documentation. 0104 dkeparam.mlprecinv = 1e-5;%Accuracy for the MatLab build-in matrix inversion. 0105 dkeparam.MUMPSparam.nproc = 1;%number of processor used by MUMPS (depends of the computer framework). Default is 1 0106 dkeparam.MUMPSparam.thresholdpivoting = '';%Relative threshold for numerical pivoting in the MUMPS partial LU decomposition (default MUMPS = 1e-2: NaN or empty or field not existing). Increasing this parameter will reduce sparcity. Less than 1. 0107 dkeparam.MUMPSparam.precinv = '';%Accuracy for the MUMPS matrix inversion.(default MUMPS: NaN or empty or field not existing). Only for the solve phase. 0108 dkeparam.PETScparam.nproc = 1;%number of processor used by PETSc (depends of the computer framework). Default is 1 0109 dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation 0110 dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs' for example. See PETSc documentation 0111 dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE. See PETSc documentation 0112 dkeparam.PETScparam.cloop = 1;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0) 0113 % 0114 dkeparam.clustermode.coll_dke_jd.scheduler.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0115 dkeparam.clustermode.coll_dke_jd.scheduler.memory = 500;%required memory (in mb) 0116 dkeparam.clustermode.eecoll_dke_yp.scheduler.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0117 dkeparam.clustermode.eecoll_dke_yp.scheduler.memory = 500;%required memory (in mb) 0118 dkeparam.clustermode.eecoll_dke_yp.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0119 dkeparam.clustermode.wave_solver_yp.scheduler.memory = 500;%required memory (in mb) 0120 % 0121 dkeparam.store_mode = 0;%Temporary storing option for sparse matrices calculations (in RAM: 0, in disk file 1) when memory space is tight 0122 dkeparam.opt_load = 0;%Load intermediate calculations 0123 dkeparam.opt_save = 0;%Save intermediate calculations 0124 % 0125 elseif strcmp(par.luke_mode,'ecrh'),%For the Electron Cyclotron wave simulations 0126 % 0127 %Mesh size calculation 0128 % 0129 dkeparam.id = 'cronos_EC'; 0130 % 0131 dkeparam.grid_uniformity = 0;%Pitch-angle and momentum grids uniformity: (0) no, (1) yes (-1 for specific test for DKE, transient and pitch-angle grid only) 0132 dkeparam.grid_method = 1;%Pitch-angle grid calculation technique: (0): old method, (1) new method 0133 dkeparam.nmhu_S = 101;%Number of angular values for the flux grid (WARNING: odd number) (>100 recommended for accurate run) 0134 dkeparam.np_S = 121;%Number of momentum values for the flux grid (pnmax/np < 0.1 highly recommended for accurate run) WARNING: must be adjusted so that effective momentum point number is larger than thepitch-angle one. 0135 dkeparam.pnmax_S = 15;%Maximum value of the flux momentum grid (pth_ref) 0136 dkeparam.rho_S = 20;% number of points for radial grid calculation 0137 dkeparam.rho0 = 0.01;%option for extra grid point near 0 (value above 0, none if 0) 0138 % 0139 %Pitch-angle grid 0140 % 0141 dkeparam.sfac0 = 1; % reduction factor for fine grid near mhu = 0 0142 dkeparam.sfac1 = 1; % reduction factor for fine grid near abs(mhu) = 1 0143 dkeparam.sfact = 1; % reduction factor for fine grid near abs(mhu) = mhut 0144 % 0145 dkeparam.snmhu0 = 5; % number of fine grid points near mhu = 0 0146 dkeparam.snmhu1 = 1; % number of fine grid points near abs(mhu) = 1 0147 dkeparam.snmhut = 1; % number of fine grid points near abs(mhu) = mhut 0148 % 0149 dkeparam.grid_power = 1.5; % power for mhu grid step varations (grid_power > 1) (Note: for grid_method = 1 only) 0150 % 0151 %Momentum grid 0152 % 0153 dkeparam.np_S_ref = 20;%position in the mpomentum grid where dpn_S is twice the minimum value(must be ajusted with snp parameter) 0154 dkeparam.snp = 1;%pgridfactor (Uniform grid: pgridfactor = 1, non-uniform grid: pgridfactor > 1) 0155 dkeparam.spfac = 20;%reduction factor for the momentum supergrid used in collision integrals 0156 % 0157 %Normalized radii at which the distribution is calculated (DKE or FP) 0158 %The syntax rho='all' may be also used if all mhu values are considered. 0159 % 0160 dkeparam.dpsin_dke = 0.01;%maximum difference between two consecutive positions used for radial derivatives 0161 % 0162 dkeparam.ref_mode = 0;%Reference values: (0) from core Te and ne, (1) from Te and ne at the normalized radius (paramater only active for a single radial value, otherwise ref_mode = 0) 0163 % 0164 % Dql and collision operator 0165 % 0166 dkeparam.coll_mode = 2;%Relativistic collision operator: (0) Relativistic Maxwellian background, (1) High-velocity limit, (2) Linearized Belaiev-Budker (3) Non-relativistic Lorentz model, (4) Non-relativistic Maxwellian background 0167 % 0168 % Ray-tracing and RF waves 0169 % 0170 dkeparam.rt_mode = 1;%RT-DKE power consistency 0171 dkeparam.abs_lim = 0;%minimum required absorption fraction 0172 dkeparam.nk_abs = 1;%maximum number of iterations on a_sdNpar for full absorption 0173 dkeparam.nit_rf = 40;%Maximum number of iteration in ray-tracing. Works only if rt_mode = 1 and more than one radial point (default = 1) 0174 dkeparam.prec_rf = 1e-2;%relative precision in RT-DKE iterations 0175 %dkeparam.n_rf_list = 1:3;%List of harmonics in the Dql calculations 0176 dkeparam.delta_rt = 0.25;%smoothering parameter of quasilinear iterations 0177 dkeparam.Dparparmax = 20;%Integration time step (1/nhuth_ref) (dtn > 0: fully implicit time differencing scheme , dtn < 0: Crank-Nicholson time differencing scheme) (WARNING: dtn less than 1 when fast electron radial transport is considered) 0178 dkeparam.npresmin = 1;%minimum number of resonant point per ray 0179 % 0180 % Run mode 0181 % 0182 dkeparam.bounce_mode = 1;%Bounce averaging option: (0) no bounce averaging, (1) bounce averaging full implicit method with 15 diagonals 0183 dkeparam.dke_mode = 0;%Fokker-Planck equation (0), Drift Kinetic equation (1) 0184 dkeparam.equil_mode = 1;%Equilibrium option: (0) circular concentric flux-surfaces, (1) numerical equilibrium 0185 dkeparam.syn_mode = 0;%Synergy calculation: (0) none, (1) yes: the code runs twice, including once without RF or electric fields. (only if dke_mode = 1) 0186 % 0187 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 0188 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) 0189 dkeparam.prec0_f = 1e-11;%Convergence level of the residu 0190 dkeparam.nit_f = 500;%Maximum number of iterations for f (> 1) 0191 dkeparam.initguess_f = 0;%Initial guess for the distribution function: (0) Maxwellian, (1) Fuchs model with Tperp for LH current drive only, (2) from a file 0192 dkeparam.initfile_f = '';%If initguess_f = 2. The file must be in the directory Project_DKE, and created by dke_4_1yp.m 0193 % 0194 dkeparam.boundary_mode_g = 0;%Number of points where the first distribution function for the Lorentz model 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) 0195 dkeparam.prec0_g = 1e-11;%Convergence level of the residu (if NaN, iterations stop with nit_g) 0196 dkeparam.nit_g = 100;%Maximum number of iterations for g 0197 % 0198 dkeparam.dtn = NaN;%Integration time step (1/nhuth_ref) (dtn > 0: fully implicit time differencing scheme , dtn < 0: Crank-Nicholson time differencing scheme) (WARNING: dtn less than 1 when fast electron radial transport is considered) 0199 % 0200 dkeparam.invproc = -1;%Iterative matrix inversion methods 0201 % (3) Matlab built-in qmr iterative method (MatLab built-in partial LU factorization used) 0202 % (2) Matlab built-in bicg iterative method (MatLab built-in partial LU factorization used) 0203 % (1) Matlab built-in cgs iterative method (MatLab built-in partial LU factorizsation used) 0204 % (0) Matlab built-in direct inversion without any LU factorization (NOT RECOMMANDED !) 0205 % (-1): MUMPS parallel (or sequential) fortran 95 direct inversion (LU matrix factorization done for each time a matrix inversion is called) 0206 % (-2): MUMPSMEX sequential fortran 95 inversion using MUMPS LU matrix factorization 0207 % (-3): SUPERLUMEX sequential direct matrix inversion (LU matrix factorization done for each time a matrix inversion is called) 0208 % (-30) PETSc interface (see data structure below) 0209 dkeparam.ludroptol = 1e-4;%Drop tolerance for the built-in MatLab partial LU decomposition (recipes: droptol = max(1/dtn/nr,param_inv(1))) 0210 dkeparam.mlmaxitinv = 4;%Maximum number of iterations for iterative built-in MatLab matrix inversion method (between 1 and 20 usualy, if empty, mlmaxitinv = min(matrixdim,20). See matlab documentation. 0211 dkeparam.mlprecinv = 1e-4;%Accuracy for the MatLab build-in matrix inversion. 0212 dkeparam.MUMPSparam.nproc = 1;%number of processor used by MUMPS (depends of the computer framework). Default is 1 0213 dkeparam.PETScparam.nproc = 3;%number of processor used by PETSc (depends of the computer framework). Default is 1 0214 dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation 0215 dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs' for example. See PETSc documentation 0216 dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE. See PETSc documentation 0217 dkeparam.PETScparam.cloop = 1;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0) 0218 % 0219 dkeparam.clustermode.coll_dke_jd.scheduler.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0220 dkeparam.clustermode.coll_dke_jd.scheduler.memory = 500;%required memory (in mb) 0221 dkeparam.clustermode.eecoll_dke_yp.scheduler.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0222 dkeparam.clustermode.eecoll_dke_yp.scheduler.memory = 500;%required memory (in mb) 0223 dkeparam.clustermode.wave_solver_yp.scheduler.mode = 0;%MatLab distributed computing environment disabled (0), enabled with the dedicated toolbox (1), enabled with a private method (2) for the function coll_dke_jd.m (MDC toolbox must be installed for option 1) 0224 dkeparam.clustermode.wave_solver_yp.scheduler.memory = 500;%required memory (in mb) 0225 % 0226 dkeparam.store_mode = 0;%Temporary storing option for sparse matrices calculations (in RAM: 0, in disk file 1) when memory space is tight 0227 dkeparam.opt_load = 0;%Load intermediate calculations 0228 dkeparam.opt_save = 0;%Save intermediate calculations 0229 % 0230 else 0231 error('The dkeparam configuration name does not exists') 0232 end 0233 % 0234 % --------- Specification par utilisateur CRONOS ------------ 0235 % 0236 % LUKE parameters 0237 % 0238 dkeparam.rt_mode = par.rt_mode;%RT-DKE power consistency 0239 dkeparam.abs_lim = par.abs_lim;%Relative level of RF power absorption that is accpeted for in the quasilinear convergence loop (0: no threshold, 1: full absorption) 0240 dkeparam.nit_rf = par.nit_rf;%Maximum number of iteration in ray-tracing. Works only if rt_mode = 1 and more than one radial point (default = 1) 0241 %dkeparam.dke_mode = par.dke_mode;%Fokker-Planck equation (0), Drift Kinetic equation (1) 0242 if isfield(par,'tn') 0243 dkeparam.tn = par.tn;% LUKE output time 0244 else 0245 dkeparam.tn = 100000; 0246 end 0247 dkeparam.invproc = par.invproc;%Iterative matrix inversion method: (1) cgs iterative method, (2) bicg iterative method, (3) qmr iterative method 0248 % 0249 if ~isfield(par,'nk_abs'), 0250 par.nk_abs = 1; 0251 end 0252 % 0253 dkeparam.nk_abs = par.nk_abs; 0254 % 0255 if ~isfield(par,'Dql') 0256 disp('Old simulation : par.Dql is missing. It is put to 20') 0257 else 0258 dkeparam.Dparparmax = par.Dql; 0259 end 0260 % 0261 % RADIAL GRID 0262 % 0263 if ischar(par.rho_S) 0264 if strcmp(par.lob,'No') 0265 eval(['par.rho_S = [',par.rho_S,'];']); 0266 else 0267 eval(['par.rho_S = [',par.rho_S,'];']); 0268 if isscalar(par.rho_S), 0269 par.rho_S=min(25,par.rho_S+5); 0270 disp(['Increase the number of radial point due to the LH negative lob : ',int2str(par.rho_S)]) 0271 end 0272 end 0273 end 0274 % 0275 if isfield(par,'rho_S') 0276 dkeparam.rho_S = par.rho_S; 0277 end 0278 % 0279 dkeparam.psin_S = []; 0280 % 0281