rundke

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 p_opt = 2;
0012 %
0013 permission = test_permissions_yp;
0014 %
0015 if ~permission 
0016     disp('Please move the script to a local folder where you have write permission before to run it')
0017     return;
0018 end
0019 %
0020 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0021 %
0022 id_simul = 'Ohm_bounce';%Simulation ID
0023 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0024 %
0025 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0026 rho_S = [linspace(0,0.2,11),linspace(0.25,0.9,14),linspace(0.925,1,5)];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0027 %
0028 id_path = '';%For all paths used by DKE solver
0029 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0030 %
0031 id_equil_cyl = 'TScyl';%For plasma equilibrium
0032 id_equil = 'TScirc_e1';%For plasma equilibrium
0033 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0036 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_display = 'NO_DISPLAY';%For output code display
0039 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 id_ohm = '';%For Ohmic electric contribution
0042 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0045 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_transpfaste = '';%For fast electron radial transport
0048 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 id_ripple = '';%For fast electron magnetic ripple losses
0051 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0052 %
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %************************************************************************************************************************************
0056 %
0057 [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);
0058 [equil_cyl] = load_structures_yp('equil',id_equil_cyl,path_equil);
0059 %
0060 %************************************************************************************************************************************
0061 %
0062 if exist('dmumpsmex');dkeparam.invproc = -2;end
0063 %
0064 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)
0065 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
0066 dkeparam.tn = [50000,100000];% 2 time steps to tn=100000 best for asymptotic solution
0067 %
0068 dkeparam.psin_S = psin_S;
0069 dkeparam.rho_S = rho_S;
0070 %
0071 epsi = 0.001;
0072 betath = 0.01;%validated for NR limit
0073 %
0074 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0075 %
0076 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0077 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0078 %equil.pzTi = 1e-10*ones(size(equil.pzTi));
0079 %
0080 equil_cyl.pTe = betath^2*mc2*ones(size(equil_cyl.pTe));
0081 equil_cyl.pzTi = betath^2*mc2*ones(size(equil_cyl.pzTi));
0082 %equil_cyl.pzTi = 1e-10*ones(size(equil_cyl.pzTi));
0083 %
0084 ohm_cyl = ohm_flat(equil_cyl,epsi);
0085 ohm = ohm_flat(equil,epsi);
0086 %
0087 dkeparam.bounce_mode = 1;
0088 %
0089 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0090 if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0091     j = Zcurr.x_tot_vfsav;
0092 else
0093     j = NaN;
0094 end  
0095 %
0096 dkeparam.bounce_mode = 0;
0097 dkeparam.mhu_S = momentumDKE.mhu_S;
0098 %
0099 [Znorm_cyl,Zcurr_cyl,ZP0_cyl,dke_out_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm_cyl,waves,transpfaste,ripple,[],[]);
0100 if dke_out_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0101     j_0 = Zcurr_cyl.x_0;
0102 else
0103     j_0 = NaN;
0104 end        
0105 %
0106 %
0107 %************************************************************************************************************************************
0108 %
0109 iar_luke = equilDKE.xrho*equilDKE.ap/equilDKE.Rp;%Circular inverse aspect ratio
0110 sigman_luke = j./j_0;
0111 %
0112 [iar,sigman_sc] = conductivity_model_yp(equil,1);%Sigmar-Coppi model (low inverse aspect ratio limit)
0113 [iar,sigman_cghk] = conductivity_model_yp(equil,2);%Connor et al. model (arbitrary inverse aspect ratio)
0114 %
0115 %************************************************************************************************************************************
0116 %
0117 figure(1),clf
0118 %
0119 leg = {'LUKE','Sigmar-Coppi','Connor et al.','Location','NorthWest'};
0120 xlim = [0,1];
0121 ylim = [1,1000];
0122 xlab = 'r/R_p';
0123 ylab = '\sigma_{cyl}/\sigma';
0124 tit = '';
0125 siz = 20+14i;
0126 %
0127 graph1D_jd(iar_luke,1./sigman_luke,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0128 graph1D_jd(iar,1./sigman_sc,0,0,'','','',NaN,xlim,ylim,'-','none','b',2,siz,gca);
0129 graph1D_jd(iar,1./sigman_cghk,0,0,'','','',leg,xlim,ylim,'-','none','g',2,siz,gca);
0130 %
0131 set(gca,'xtick',[0:0.2:1])
0132 set(gca,'ytick',[1 10 100 1000 10000])
0133 set(gca,'YMinorGrid','off')
0134 set(gca,'YMinorTick','on')
0135 %
0136 print_jd(p_opt,'fig_sigma_bounce','./figures',1)
0137 %
0138 %************************************************************************************************************************************
0139 %
0140 %
0141 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0142 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);