rundke_val

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 = 'Ohm_bounce_val';%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 = [];%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_cyl = 'TScyl';%For plasma equilibrium
0030 id_equil = 'TScirc_e1';%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 = 'NONUNIFORM10010020';%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 = 'PARTIAL_VISUAL';%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 [equil_cyl] = load_structures_yp('equil',id_equil_cyl,path_equil);
0057 %
0058 %************************************************************************************************************************************
0059 %
0060 if exist('dmumpsmex');dkeparam.invproc = -2;end
0061 %
0062 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)
0063 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
0064 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0065 %
0066 dkeparam.psin_S = psin_S;
0067 %
0068 iar_list = [0.1,0.2,0.3,0.5,0.8];
0069 %
0070 epsi = 0.001;
0071 betath = 0.01;%validated for NR limit
0072 %
0073 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0074 %
0075 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0076 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0077 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0078 %
0079 equil_cyl.pTe = betath^2*mc2*ones(size(equil_cyl.pTe));
0080 equil_cyl.pzTi = betath^2*mc2*ones(size(equil_cyl.pzTi));
0081 %equil_cyl.pzTi = 1e-10*ones(size(equil_cyl.pzTi));
0082 %
0083 ohm_cyl = ohm_flat(equil_cyl,epsi);
0084 ohm = ohm_flat(equil,epsi);
0085 %
0086 ap = equil.ptx(end,1);
0087 rho_list = iar_list*equil.Rp/ap;
0088 %
0089 nrho = length(rho_list);
0090 j_0 = NaN(1,nrho);
0091 j = NaN(1,nrho);
0092 %
0093 for irho = 1:nrho,
0094     %
0095     dkeparam.bounce_mode = 1;
0096     dkeparam.mhu_S = [];
0097     %
0098     dkeparam.rho_S = rho_list(irho);
0099     %
0100     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0101     %
0102     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0103         j(irho) = Zcurr.x_tot_vfsav;
0104     end  
0105     %
0106     dkeparam.bounce_mode = 0;
0107     dkeparam.mhu_S = momentumDKE.mhu_S;
0108     %
0109     dkeparam.rho_S = 0.5;
0110     %
0111     [Znorm_cyl,Zcurr_cyl,ZP0_cyl,dke_out_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm_cyl,waves,transpfaste,ripple,[],[]);
0112     %
0113     if dke_out_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0114         j_0(irho) = Zcurr_cyl.x_0;
0115     end  
0116     %
0117 end
0118 %
0119 %************************************************************************************************************************************
0120 %
0121 sigman_luke = j./j_0;
0122 %
0123 [iar,sigman_sc] = conductivity_model_yp(equil,1,rho_list);%Sigmar-Coppi model (low inverse aspect ratio limit)
0124 [iar,sigman_cghk] = conductivity_model_yp(equil,2,rho_list);%Connor et al. model (arbitrary inverse aspect ratio)
0125 %
0126 %************************************************************************************************************************************
0127 %
0128 format
0129 %
0130 delete res_ohm_val
0131 %
0132 diary res_ohm_val
0133 %
0134 disp(['Comparison LUKE/Coppi/Connor - NR calculations (betath = 0.01)'])
0135 disp(['----------------------'])
0136 disp(['--> eps = 0.1 : sigma/sigma_cyl = ',num2str(sigman_luke(1)),'/',num2str(sigman_sc(1)),'/',num2str(sigman_cghk(1))]) 
0137 disp(['--> eps = 0.2 : sigma/sigma_cyl = ',num2str(sigman_luke(2)),'/',num2str(sigman_sc(2)),'/',num2str(sigman_cghk(2))]) 
0138 disp(['--> eps = 0.3 : sigma/sigma_cyl = ',num2str(sigman_luke(3)),'/',num2str(sigman_sc(3)),'/',num2str(sigman_cghk(3))]) 
0139 disp(['--> eps = 0.5 : sigma/sigma_cyl = ',num2str(sigman_luke(4)),'/',num2str(sigman_sc(4)),'/',num2str(sigman_cghk(4))]) 
0140 disp(['--> eps = 0.8 : sigma/sigma_cyl = ',num2str(sigman_luke(5)),'/',num2str(sigman_sc(5)),'/',num2str(sigman_cghk(5))]) 
0141 %
0142 diary off
0143 %
0144 %************************************************************************************************************************************
0145 %
0146 %
0147 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0148 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);