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 opt_proc = 0;%test for numerical performance test with number of processors
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 = 'LH_karney_transp_bounce_time';%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:0.025:0.65,0.7:0.1:1];%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 = 'JETh77';%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 = 'NONUNIFORM10010020';%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 = 'PARTIAL_VISUAL';%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 = 'magneticturbulence';%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 wavestruct.omega_lh = [4]*2*pi*1e9; %(GHz -> rad/s). Wave frequency [1,1] Indicative, no effect in small FLR limit opt_lh > 0
0061 %Option parameter for cross-comparison between old LH code:
0062 %    - (1): 1/vpar dependence
0063 %    - (2): no 1/vpar dependence and old grid technique for Dlh calculations (Karney, Shoucri) (see rfdiff_dke_jd)
0064 wavestruct.opt_lh = 2; % [1,1]
0065 %
0066 % Choose (vparmin_lh,vparmax_lh) or (Nparmin_lh,Nparmax_lh) for square n// LH wave power spectrum,
0067 % or (Npar_lh,dNpar_lh) for Gaussian shape
0068 %
0069 wavestruct.norm_ref = 1;%Normalization procedure for the LH quasilinear diffusion coefficient and spectrum boundaries
0070 %
0071 wavestruct.yvparmin_lh = [4];%LH wave square N// Spectrum: Lower limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0072 wavestruct.yvparmax_lh = [7];%LH wave square N// Spectrum: Upper limit of the plateau (vth_ref or vth) [1,n_scenario_lh]
0073 %
0074 wavestruct.yNparmin_lh = [NaN];%LH wave square N// Spectrum: Lower limit [1,n_scenario_lh]
0075 wavestruct.yNparmax_lh = [NaN];%LH wave square N// Spectrum: Upper limit [1,n_scenario_lh]
0076 wavestruct.yNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: peak [1,n_scenario_lh]
0077 wavestruct.ydNpar_lh = [NaN];%LH wave Gaussian N// Spectrum: width [1,n_scenario_lh]
0078 %
0079 %   Note: this diffusion coefficient is different from the general QL D0. It has a benchmarking purpose only
0080 wavestruct.yD0_in_c_lh = [2];%Central LH QL diffusion coefficient (nhuth_ref*pth_ref^2 or nhuth*pth^2) [1,n_scenario_lh]
0081 %
0082 wavestruct.yD0_in_lh_prof = [1];%Quasilinear diffusion coefficient radial profile: (0) uniform, (1) gaussian radial profile [1,n_scenario_lh]
0083 wavestruct.ypeak_lh = [0.4];%Radial peak position of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0084 wavestruct.ywidth_lh = [0.1];%Radial width of the LH quasi-linear diffusion coefficient (r/a on midplane) [1,n_scenario_lh]
0085 %
0086 wavestruct.ythetab_lh = [0]*pi/180;%(deg -> rad). Poloidal location of LH beam [0..2pi] [1,n_scenario_lh]
0087 %               (0) from local values Te and ne, (1) from central values Te0 and ne0
0088 %
0089 %************************************************************************************************************************************
0090 %
0091 if exist('dmumpsmex');dkeparam.invproc = -2;end
0092 %
0093 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)
0094 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
0095 dkeparam.dtn = 10;%Normalized integration time step
0096 dkeparam.nit_f = 1;%Maximum number of internal iterations for f (>= 1)
0097 %
0098 dkeparam.psin_S = psin_S;
0099 dkeparam.rho_S = rho_S;
0100 %
0101 ntime_step = 10;%number of time steps
0102 %
0103 transpfaste.Dr0 = 0;
0104 %
0105 % Initial magnetic equilibrium (JETh77)
0106 %
0107 R_sep(1) = 3.0;%  Geometrical major radius of the separatrix (m)
0108 a_sep(1) = 0.9;%  Geometrical minor radius of the separatrix (m)
0109 Z_sep(1) = 0.0;%  Geometrical vertical shift of the separatrix (m)
0110 %
0111 xpoint(1,1) = 1.687;%Upper altitude X point in minor radius unit
0112 xpoint(2,1) = 0.466;%Upper triangularity in minor radius unit
0113 xpoint(3,1) = 0;%Upper separatrix angle (R,X)  (LFS, degrees)
0114 xpoint(4,1) = 0;%Upper separatrix angle (-R,X) (HFS, degrees)
0115 xpoint(5,1) = 2.001;%Lower altitude X point in minor radius unit (always positive)
0116 xpoint(6,1) = 0.568;%Lower triangularity in minor radius unit
0117 xpoint(7,1) = 22.46;%Lower separatrix angle (R,X)  (LFS, degrees)
0118 xpoint(8,1) = 67.92;%Lower separatrix angle (-R,X)  (HFS, degrees)
0119 %
0120 npsi(1) = 101;%ONLY 101 radial profile grid points for HELMEX77
0121 ntheta(1) = 65;%ONLY 65 poloidal grid points for HELMEX77
0122 %
0123 Ip(1) = 2;%  Plasma current (MA)
0124 %
0125 qmin(1) = 1;%  Safety factor q0 at plasma center
0126 eq(1) = 3;%  Exponent for q radial profile (q = (qmax - qmin)*(r/a).^eq + qmin, qmax calaculated by the Ampere's theorem at the plasma edge with a circular plasma cross-section)
0127 qopt(1) = 1;%Option for q profile. (0): constant (default, uniform current, psi=rho^2), (1): from qmin and eq
0128 %
0129 Bt(1) = 3.5;%  Toroidal magnetic field on the magnetic axis (T)
0130 %
0131 % Kinetic parameters (profiles)
0132 %
0133 Zi(1,:) =  [1,1,1,6];%  Ion types: (1) H/D/T, (2) He, ..., (6) C [1,p] (WARNING: Zi must be [1,1,1,imp1,imp2] for hydrogen plasmas)
0134 mi(1,:) =  [1,2,3,12];%  Ion mass (uma) [1,p] (WARNING: Zi must be [1,2,3,mimp1,mimp2] for hydrogen plasmas)
0135 fi(1,:) =  [0,1,0];%  Hydrogen isotopic fraction (H/D/T) [1,3] (WARNING: only used when hydrogen plasmas are considered)
0136 %
0137 Te0(1) = 5;%  Core electron temperature (keV)
0138 Tea(1) = 0.5;%  Edge electron temperature (keV)
0139 eTe1(1) = 2;%  Exponent #1 for Te profile
0140 eTe2(1) = 2;%  Exponent #2 for Te profile
0141 %
0142 ne0(1) = 2.0e19;%  Core electron density (m-3)
0143 nea(1) = 2.0e17;%  Edge electron density (m-3)
0144 ene1(1) = 2;%  Exponent #1 for ne profile
0145 ene2(1) = 0.5;%  Exponent #2 for ne profile
0146 %
0147 Ti0(1) = 2;%  Core ion temperature (keV)
0148 Tia(1) = 0.2;%  Edge ion temperature (keV)
0149 eTi1(1) = 2;%  Exponent #1 for Ti profile
0150 eTi2(1) = 2;%  Exponent #2 for Ti profile
0151 %
0152 Zeff0(1) = 1.5;%  Core effective charge (a.u.)
0153 Zeffa(1) = 1.5;%  Edge effective charge (a.u.)
0154 eZeff(1) = 0;%  Define the slope at the radial point where the derivative is maximum (must be >> 1)
0155 xZeff(1) = 0.5;%Radial normalized position at which the derivative is maximum
0156 %
0157 [ppsin,ppsi_apRp,theta,x,y,Bx,By,BPHI,pBpp,pq_Rpap,pj,pmag,Ip_test] = idealequilcyl_yp(a_sep(1),R_sep(1),Z_sep(1),Bt(1),Ip(1),qmin(1)*R_sep(1)/a_sep(1),eq(1),qopt(1),npsi(1),ntheta(1));%Cylindrical magnetic equilibrium
0158 %
0159 [prho,pTe,pne,pzTi,pzni,zZi,zmi,fi,pkin] = idealprof_yp(Zi(1,:),mi(1,:),fi(1,:),Te0(1),Tea(1),eTe1(1),eTe2(1),ne0(1),nea(1),ene1(1),ene2(1),Ti0(1),Tia(1),eTi1(1),eTi2(1),Zeff0(1),Zeffa(1),eZeff(1),xZeff(1),npsi(1));%Profiles
0160 %
0161 pj_min = min(pj); 
0162 if pj_min < 0,pj = pj - pj_min;end,%pj must be non-negative. Only the shape is considered and Ip for the total plasma current
0163 [equil_magnetic,prho_out,pspsin_out,pspsinT_out,pp_out] = equil_magnetic_helmex77_yp(a_sep(1),R_sep(1),Z_sep(1),xpoint(:,1),Bt(1),Ip(1),ppsi_apRp*R_sep(1)/a_sep(1),pj,pkin,ntheta(1),0);
0164 %
0165 equil_prof.pTe = pTe;
0166 equil_prof.pne = pne;
0167 equil_prof.pzTi = pzTi;
0168 equil_prof.pzni = pzni;
0169 equil_prof.zZi = zZi;
0170 equil_prof.zmi = zmi;
0171 %
0172 clear equil
0173 %
0174 equil{1} = conc_struct_jd(equil_magnetic,equil_prof);
0175 equil{1}.id = id_equil;
0176 %
0177 equilibrium_jd(equil{1},'',2,'linear',50);
0178 drawnow
0179 %
0180 % Magnetic equilibrium independent of time (fixed equilibrium)
0181 %
0182 waves{1} = make_idealLHwave_jd(equil{1},wavestruct);
0183 %
0184 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul,dkepath,equil{1},dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0185 %
0186 j_fe(:,1) = Zcurr.x_0;
0187 P_fe(:,1) = ZP0.x_rf_fsav;
0188 Itot_fe(1) = sum(Zcurr.x_0.*equilDKE.xdA_dke*mksa.j_ref);
0189 Ptot_2piRp_fe(1) = sum(ZP0.x_rf_fsav.'.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0190 %
0191 for it = 2:ntime_step,
0192     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul,dkepath,equil{1},dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out.Zf0_interp,[]);
0193     %
0194     j_fe(:,it) = Zcurr.x_0;
0195     P_fe(:,it) = ZP0.x_rf_fsav;
0196     Itot_fe(it) = sum(Zcurr.x_0.*equilDKE.xdA_dke*mksa.j_ref);
0197     Ptot_2piRp_fe(it) = sum(ZP0.x_rf_fsav.'.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0198 end
0199 %
0200 % Magnetic equilibrium evolves with time (variable equilibrium)
0201 %
0202 waves{1} = make_idealLHwave_jd(equil{1},wavestruct);
0203 %
0204 [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul,dkepath,equil{1},dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0205 %
0206 j_ve(:,1) = Zcurr.x_0;
0207 P_ve(:,1) = ZP0.x_rf_fsav;
0208 Itot_ve(1) = sum(Zcurr.x_0.*equilDKE.xdA_dke*mksa.j_ref);
0209 Ptot_2piRp_ve(1) = sum(ZP0.x_rf_fsav.'.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0210 %
0211 % Other time steps
0212 %
0213 for it = 2:ntime_step,
0214     R_sep(it) = 3.0 - 1.5*(it/ntime_step);%  Geometrical major radius of the separatrix (m)
0215     a_sep(it) = 0.9;%  Geometrical minor radius of the separatrix (m)
0216     Z_sep(it) = 0.0;%  Geometrical vertical shift of the separatrix (m)
0217     %
0218     xpoint(1,it) = 1.687;%Upper altitude X point in minor radius unit
0219     xpoint(2,it) = 0.466;%Upper triangularity in minor radius unit
0220     xpoint(3,it) = 0;%Upper separatrix angle (R,X)  (LFS, degrees)
0221     xpoint(4,it) = 0;%Upper separatrix angle (-R,X) (HFS, degrees)
0222     xpoint(5,it) = 2.001;%Lower altitude X point in minor radius unit (always positive)
0223     xpoint(6,it) = 0.568 - 0.102*(it/ntime_step);%Lower triangularity in minor radius unit
0224     xpoint(7,it) = 22.46 - 22.46*(it/ntime_step);%Lower separatrix angle (R,X)  (LFS, degrees)
0225     xpoint(8,it) = 67.92 - 67.92*(it/ntime_step);%Lower separatrix angle (-R,X)  (HFS, degrees)
0226     %
0227     npsi(it) = 101;%ONLY 101 radial profile grid points for HELMEX77
0228     ntheta(it) = 65;%ONLY 65 poloidal grid points for HELMEX77
0229     %
0230     Ip(it) = 2;%  Plasma current (MA)
0231     %
0232     qmin(it) = 1;%  Safety factor q0 at plasma center
0233     eq(it) = 3;%  Exponent for q radial profile (q = (qmax - qmin)*(r/a).^eq + qmin, qmax calaculated by the Ampere's theorem at the plasma edge with a circular plasma cross-section)
0234     qopt(it) = 1;%Option for q profile. (0): constant (default, uniform current, psi=rho^2), (1): from qmin and eq
0235     %
0236     Bt(it) = 3.5;%  Toroidal magnetic field on the magnetic axis (T)
0237     %
0238     % Kinetic parameters (profiles)
0239     %
0240     Zi(it,:) =  [1,1,1,6];%  Ion types: (1) H/D/T, (2) He, ..., (6) C [1,p] (WARNING: Zi must be [1,1,1,imp1,imp2] for hydrogen plasmas)
0241     mi(it,:) =  [1,2,3,12];%  Ion mass (uma) [1,p] (WARNING: Zi must be [1,2,3,mimp1,mimp2] for hydrogen plasmas)
0242     fi(it,:) =  [0,1,0];%  Hydrogen isotopic fraction (H/D/T) [1,3] (WARNING: only used when hydrogen plasmas are considered)
0243     %
0244     Te0(it) = 5;%  Core electron temperature (keV)
0245     Tea(it) = 0.5;%  Edge electron temperature (keV)
0246     eTe1(it) = 2;%  Exponent #1 for Te profile
0247     eTe2(it) = 2;%  Exponent #2 for Te profile
0248     %
0249     ne0(it) = 2.0e19;%  Core electron density (m-3)
0250     nea(it) = 2.0e17;%  Edge electron density (m-3)
0251     ene1(it) = 2;%  Exponent #1 for ne profile
0252     ene2(it) = 0.5;%  Exponent #2 for ne profile
0253     %
0254     Ti0(it) = 2;%  Core ion temperature (keV)
0255     Tia(it) = 0.2;%  Edge ion temperature (keV)
0256     eTi1(it) = 2;%  Exponent #1 for Ti profile
0257     eTi2(it) = 2;%  Exponent #2 for Ti profile
0258     %
0259     Zeff0(it) = 1.5;%  Core effective charge (a.u.)
0260     Zeffa(it) = 1.5;%  Edge effective charge (a.u.)
0261     eZeff(it) = 0;%  Define the slope at the radial point where the derivative is maximum (must be >> 1)
0262     xZeff(it) = 0.5;%Radial normalized position at which the derivative is maximum
0263     %
0264     [ppsin,ppsi_apRp,theta,x,y,Bx,By,BPHI,pBpp,pq_Rpap,pj,pmag,Ip_test] = idealequilcyl_yp(a_sep(it),R_sep(it),Z_sep(it),Bt(it),Ip(it),qmin(it)*R_sep(it)/a_sep(it),eq(it),qopt(it),npsi(it),ntheta(it));%Cylindrical magnetic equilibrium
0265     %
0266     [prho,pTe,pne,pzTi,pzni,zZi,zmi,fi,pkin] = idealprof_yp(Zi(it,:),mi(it,:),fi(it,:),Te0(it),Tea(it),eTe1(it),eTe2(it),ne0(it),nea(it),ene1(it),ene2(it),Ti0(it),Tia(it),eTi1(it),eTi2(it),Zeff0(it),Zeffa(it),eZeff(it),xZeff(it),npsi(it));%Profiles
0267     %
0268     pj_min = min(pj); 
0269     if pj_min < 0,pj = pj - pj_min;end,%pj must be non-negative. Only the shape is considered and Ip for the total plasma current
0270     [equil_magnetic,prho_out,pspsin_out,pspsinT_out,pp_out] = equil_magnetic_helmex77_yp(a_sep(it),R_sep(it),Z_sep(it),xpoint(:,it),Bt(it),Ip(it),ppsi_apRp*R_sep(it)/a_sep(it),pj,pkin,ntheta(it),0);
0271     %
0272     equil_prof.pTe = pTe;
0273     equil_prof.pne = pne;
0274     equil_prof.pzTi = pzTi;
0275     equil_prof.pzni = pzni;
0276     equil_prof.zZi = zZi;
0277     equil_prof.zmi = zmi;
0278     %
0279     equil{it} = conc_struct_jd(equil_magnetic,equil_prof);
0280     equil{it}.id = id_equil;
0281     %
0282     equilibrium_jd(equil{it},'',2,'linear',50);
0283     drawnow
0284     %
0285     waves{1} = make_idealLHwave_jd(equil{it},wavestruct);
0286     %
0287     [Znorm,Zcurr,ZP0,dke_out,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa] = main_dke_yp(id_simul,dkepath,equil{it},dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,dke_out.Zf0_interp,[]);
0288     %
0289     j_ve(:,it) = Zcurr.x_0*mksa.j_ref;
0290     P_ve(:,it) = ZP0.x_rf_fsav*mksa.P_ref;
0291     Itot_ve(it) = sum(Zcurr.x_0.*equilDKE.xdA_dke*mksa.j_ref);
0292     Ptot_2piRp_ve(it) = sum(ZP0.x_rf_fsav.'.*equilDKE.xdV_2piRp_dke*mksa.P_ref);
0293 end
0294 %
0295 Eta_fe = equilDKE.xne(1)*Itot_fe./Ptot_2piRp_fe/(2*pi)/1e19;
0296 Eta_ve = equilDKE.xne(1)*Itot_ve./Ptot_2piRp_fe/(2*pi)/1e19;
0297 %
0298 strcolor = ['b','g','r','c','m','y','k'];
0299 figure(1000),
0300 for it = 1:ntime_step,
0301     graph1D_jd(equilDKE.xrho,j_fe(:,it),0,0,'\rho','Current density (MA.m-2)','Fixed plasma shape',NaN,'','','-','none',strcolor(rem(it,length(strcolor))+1),2);
0302 end
0303 %
0304 figure(1001),
0305 for it = 1:ntime_step,
0306     graph1D_jd(equilDKE.xrho,j_ve(:,it),0,0,'\rho','Current density (MA.m-2)','Time evolving plasma shape',NaN,'','','-','none',strcolor(rem(it,length(strcolor))+1),2);
0307 end
0308 %
0309 figure(1002),
0310 graph1D_jd([1:ntime_step]*dkeparam.dtn,abs(Itot_fe),0,0,'t*\nu','Current (MA)','JETh77 @ t = 0',NaN,'','','-','none','b',2);
0311 graph1D_jd([1:ntime_step]*dkeparam.dtn,abs(Itot_ve),0,0,'t*\nu','Current (MA)','JETh77 @ t = 0',NaN,[1,ntime_step]*dkeparam.dtn,'','--','none','r',2);
0312 legend('Fixed','Evolving')
0313 %
0314 figure(1003),
0315 graph1D_jd([1:ntime_step]*dkeparam.dtn,abs(Ptot_2piRp_fe),0,0,'t*\nu','Power (MW)','JETh77 @ t = 0',NaN,'','','-','none','b',2);
0316 graph1D_jd([1:ntime_step]*dkeparam.dtn,abs(Ptot_2piRp_ve),0,0,'t*\nu','Power (MW)','JETh77 @ t = 0',NaN,[1,ntime_step]*dkeparam.dtn,'','--','none','r',2);
0317 legend('Fixed','Evolving')
0318 %
0319 figure(1004),
0320 graph1D_jd([1:ntime_step]*dkeparam.dtn,Eta_fe,0,0,'t*\nu','\eta (A/W/m2)','JETh77 @ t = 0',NaN,'','','-','none','b',2);
0321 graph1D_jd([1:ntime_step]*dkeparam.dtn,Eta_ve,0,0,'t*\nu','\eta (A/W/m2)','JETh77 @ t = 0',NaN,[1,ntime_step]*dkeparam.dtn,'','--','none','r',2);
0322 legend('Fixed','Evolving')
0323 %
0324 %************************************************************************************************************************************
0325 %
0326 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0327 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);