rundke_val

PURPOSE ^

Script for running the DKE solver (can be modified by the user for specific simulations)

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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 = 'Runaway_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 = 1;%Local normalization of f0 at each iteration: (0) no, the default value when the numerical conservative scheme is correct, (1) yes
0064 dkeparam.tn = 10000;%time for asymptotic solution with norm_mode_f = 1
0065 dkeparam.dtn = 1000;%10 time steps required for accurate runaway solution - see rundke_dtn
0066 %
0067 dkeparam.psin_S = psin_S;
0068 iar_list = [0.1,0.2,0.3,0.5,0.8];
0069 %
0070 epsi = 0.04;
0071 betath = 0.001;%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 dkeparam.np_S = 201;
0087 dkeparam.nmhu_S = 201;
0088 %
0089 ap = equil.ptx(end,1);
0090 rho_list = iar_list*equil.Rp/ap;
0091 %
0092 nrho = length(rho_list);
0093 RR_0 = NaN(1,nrho);
0094 RR = NaN(1,nrho);
0095 %
0096 for irho = 1:nrho,
0097     %
0098     dkeparam.bounce_mode = 1;
0099     dkeparam.mhu_S = [];
0100     %
0101     dkeparam.rho_S = rho_list(irho);
0102     %
0103     [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,[],[]);
0104     if dke_out.residu_f{end}(end) <= dkeparam.prec0_f,
0105         RR(irho) = dke_out.XxRR_fsav(end,:);
0106     end  
0107     %
0108     dkeparam.bounce_mode = 0;
0109     dkeparam.mhu_S = momentumDKE.mhu_S;
0110     %
0111     dkeparam.rho_S = 0.5;
0112     %
0113     [Znorm_cyl,Zcurr_cyl,ZP0_cyl,dke_out_cyl] = main_dke_yp(id_simul,dkepath,equil_cyl,dkeparam,dkedisplay,ohm_cyl,waves,transpfaste,ripple,[],[]);
0114     if dke_out_cyl.residu_f{end}(end) <= dkeparam.prec0_f,
0115         RR_0(irho) = dke_out_cyl.XxRR_fsav(end,:);
0116     end        
0117 end
0118 %
0119 %************************************************************************************************************************************
0120 %
0121 RRn = RR./RR_0;
0122 %
0123 %************************************************************************************************************************************
0124 %
0125 format
0126 %
0127 delete res_runaway_val
0128 %
0129 diary res_runaway_val
0130 %
0131 disp(['Comparison LUKE/Coppi/Connor - NR calculations (betath = 0.001)'])
0132 disp(['----------------------'])
0133 disp(['--> eps = 0.1 : RR/RR_cyl = ',num2str(RRn(1))]) 
0134 disp(['--> eps = 0.2 : RR/RR_cyl = ',num2str(RRn(2))]) 
0135 disp(['--> eps = 0.3 : RR/RR_cyl = ',num2str(RRn(3))]) 
0136 disp(['--> eps = 0.5 : RR/RR_cyl = ',num2str(RRn(4))]) 
0137 disp(['--> eps = 0.8 : RR/RR_cyl = ',num2str(RRn(5))]) 
0138 %
0139 diary off
0140 %
0141 %************************************************************************************************************************************
0142 %
0143 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0144 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0145 
0146 
0147
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.