rundke_conv

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_conv';%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 = [0.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 = 'TScyl';%For plasma equilibrium
0032 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0033 %
0034 id_dkeparam = 'UNIFORM10010020';%For DKE code parameters
0035 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0036 %
0037 id_display = 'NO_DISPLAY';%For output code display
0038 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0039 %
0040 id_ohm = '';%For Ohmic electric contribution
0041 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0042 %
0043 ids_wave = {''};%For RF waves contribution (put all the type of waves needed)
0044 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0045 %
0046 id_transpfaste = '';%For fast electron radial transport
0047 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0048 %
0049 id_ripple = '';%For fast electron magnetic ripple losses
0050 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0051 %
0052 %************************************************************************************************************************************
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %
0056 [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);
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.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
0064 %
0065 dkeparam.psin_S = psin_S;
0066 dkeparam.rho_S = rho_S;
0067 %
0068 betath = 0.01;%validated for NR limit
0069 %
0070 epsi = 0.001;
0071 %
0072 [qe,me,mp,mn,e0,mu0,re,mc2] = pc_dke_yp;%Universal physics constants
0073 %
0074 equil.pTe = betath^2*mc2*ones(size(equil.pTe));
0075 equil.pzTi = betath^2*mc2*ones(size(equil.pzTi));
0076 %
0077 ohm = ohm_flat(equil,epsi);
0078 %
0079 sigma_Karney_nr_2 = 7.429;
0080 %
0081 % Building the distribution function to a quasi steady-state for the
0082 % ohmic problem : tn = 100000/dtn = 1.
0083 %
0084 nit_f = 50;
0085 %
0086 dkeparam.tn = 100000;
0087 dkeparam.dtn = NaN;%specified by tn
0088 dkeparam.nit_f = nit_f;
0089 dkeparam.prec0_f = -1;%to reach end of Legendre iterations
0090 %
0091 [dummy,dummy,dummy,dke_out] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0092 sigma_2 = dke_out.curr0{end}/epsi;
0093 %
0094 sigma_Z1s = 1/100*sqrt(me*qe/10)/(1.65*mu0*e0^2);%Spitzer conductivity (from 'Tokamaks' by J. Wesson, Oxford Science Publication, 3rd Eds., 2004)
0095 %
0096 dcurr_2 = abs(dke_out.curr0{end} - dke_out.curr0{end}(end))/dke_out.curr0{end}(end);
0097 dnorm_2 = abs(dke_out.normf0{end} - dke_out.normf0{end}(end))/dke_out.normf0{end}(end);
0098 %
0099 %************************************************************************************************************************************
0100 %
0101 figure(1),clf
0102 %
0103 leg = {'LUKE','Karney','Spitzer'};
0104 xlim = [0,nit_f];
0105 ylim = [0,10];
0106 xlab = '# iterations';
0107 ylab = '\sigma';
0108 tit = '';
0109 siz = 20+14i;
0110 %
0111 graph1D_jd(0:nit_f,sigma_2,0,0,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0112 graph1D_jd(xlim,[sigma_Karney_nr_2,sigma_Karney_nr_2],0,0,'','','',NaN,xlim,ylim,'--','none','r',2,siz,gca);
0113 graph1D_jd(xlim,[sigma_Z1s,sigma_Z1s],0,0,'','','',leg,xlim,ylim,'--','none','k',2,siz,gca);
0114 %
0115 set(gca,'ytick',[0:2:10])
0116 %
0117 figure(2),clf
0118 %
0119 ylim = 10.^[-24,-8];
0120 ylab = 'R_f';
0121 %
0122 graph1D_jd(1:nit_f,dke_out.residu_f{end},0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0123 %
0124 set(gca,'xtick',[0:0.2:1]*xlim(2))
0125 set(gca,'ytick',[1e-24 1e-20 1e-16 1e-12 1e-8])
0126 set(gca,'YMinorGrid','off')
0127 set(gca,'YMinorTick','on')
0128 %
0129 figure(3),clf
0130 %
0131 ylim = 10.^[-20,0];
0132 ylab = '(j-j_f)/j_f';
0133 %
0134 graph1D_jd(0:nit_f,dcurr_2,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0135 %
0136 set(gca,'xtick',[0:0.2:1]*xlim(2))
0137 set(gca,'ytick',[1e-20 1e-15 1e-10 1e-5 1])
0138 set(gca,'YMinorGrid','off')
0139 set(gca,'YMinorTick','on')
0140 %
0141 figure(4),clf
0142 %
0143 ylim = 10.^[-20,0];
0144 ylab = '(n-n_f)/n_f';
0145 %
0146 graph1D_jd(0:nit_f,dnorm_2,0,1,xlab,ylab,tit,NaN,xlim,ylim,'-','none','r',2,siz,gca,0.9,0.7,0.7);
0147 %
0148 set(gca,'xtick',[0:0.2:1]*xlim(2))
0149 set(gca,'ytick',[1e-20 1e-15 1e-10 1e-5 1])
0150 set(gca,'YMinorGrid','off')
0151 set(gca,'YMinorTick','on')
0152 %
0153 % print_jd(2,'fig_j_conv','./figures',1)
0154 % print_jd(2,'fig_Rf_conv','./figures',2)
0155 % print_jd(2,'fig_jn_conv','./figures',3)
0156 % print_jd(2,'fig_nn_conv','./figures',4)
0157 %
0158 print_jd(2,'fig_ohm_sigma_conv_ss','./figures',1)
0159 print_jd(2,'fig_ohm_Rf_conv_ss','./figures',2)
0160 print_jd(2,'fig_ohm_jn_conv_ss','./figures',3)
0161 print_jd(2,'fig_ohm_nn_conv_ss','./figures',4)
0162 %
0163 %************************************************************************************************************************************
0164 %
0165 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0166 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);