diff --git a/matlab/setup.m b/matlab/setup.m index 0a37604..cefc6bb 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,142 +1,143 @@ %% ________________________________________________________________________ SIMDIR = ['../results/',SIMID,'/']; % Grid parameters GRID.pmaxe = PMAXE; % Electron Hermite moments GRID.jmaxe = JMAXE; % Electron Laguerre moments GRID.pmaxi = PMAXI; % Ion Hermite moments GRID.jmaxi = JMAXI; % Ion Laguerre moments GRID.Nx = N; % x grid resolution GRID.Lx = L; % x length GRID.Ny = N * (1-KXEQ0) + KXEQ0; % y '' GRID.Ly = L * (1-KXEQ0); % y '' GRID.Nz = Nz; % z '' GRID.Lz = Lz; % z '' GRID.kpar = KPAR; % Model parameters MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CLOS = CLOS; MODEL.NL_CLOS = NL_CLOS; if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end; MODEL.mu = MU; MODEL.mu_p = MU_P; MODEL.mu_j = MU_J; MODEL.nu = NU; % hyper diffusive coefficient nu for HW % temperature ratio T_a/T_e MODEL.tau_e = TAU; MODEL.tau_i = TAU; % mass ratio sqrt(m_a/m_i) MODEL.sigma_e = 0.0233380; MODEL.sigma_i = 1.0; % charge q_a/e MODEL.q_e =-1.0; MODEL.q_i = 1.0; if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end; % gradients L_perp/L_x MODEL.eta_n = ETAN; % source term kappa for HW MODEL.eta_T = ETAT; % Temperature MODEL.eta_B = ETAB; % Magnetic MODEL.lambdaD = LAMBDAD; % if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KX0KH),'_A_',num2str(A0KH)]; end; % Time integration and intialization parameters TIME_INTEGRATION.numerical_scheme = '''RK4'''; if (INIT_PHI); INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end; INITIAL.INIT_ZF = INIT_ZF; INITIAL.init_background = (INIT_ZF>0)*ZF_AMP; INITIAL.init_noiselvl = NOISE0; INITIAL.iseed = 42; INITIAL.mat_file = '''null'''; if (abs(CO) == 2) %Sugama operator INITIAL.mat_file = ['''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''']; elseif (abs(CO) == 3) %pitch angle operator INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''']; elseif (CO == 4) % Full Coulomb GK disp('Warning, FCGK not implemented yet') elseif (CO == -1) % DGDK disp('Warning, DGDK not debugged') end % Naming and creating input file if (CO == -3); CONAME = 'PADK'; elseif(CO == -2); CONAME = 'SGDK'; elseif(CO == -1); CONAME = 'DGDK'; elseif(CO == 0); CONAME = 'LB'; elseif(CO == 1); CONAME = 'DGGK'; elseif(CO == 2); CONAME = 'SGGK'; elseif(CO == 3); CONAME = 'PAGK'; end if (CLOS == 0); CLOSNAME = 'Trunc.'; elseif(CLOS == 1); CLOSNAME = 'Clos. 1'; elseif(CLOS == 2); CLOSNAME = 'Clos. 2'; end if (PMAXE == PMAXI) && (JMAXE == JMAXI) degngrad = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)]; else degngrad = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),... '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)]; end degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',... - CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e']; + CONAME,'_mu_%0.0e']; degngrad = sprintf(degngrad,[NU,MU]); if ~NON_LIN; degngrad = ['lin_',degngrad]; end if (Nz == 1) resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'_']; + gridname = ['L_',num2str(L),'_']; else resolution = [num2str(GRID.Nx),'x',num2str(GRID.Ny/2),'x',num2str(GRID.Nz),'_']; + gridname = ['L_',num2str(L),'_Lz_',num2str(Lz),'_']; end -gridname = ['L_',num2str(L),'_']; if (exist('PREFIX','var') == 0); PREFIX = []; end; if (exist('SUFFIX','var') == 0); SUFFIX = []; end; PARAMS = [PREFIX,resolution,gridname,degngrad,SUFFIX]; BASIC.RESDIR = [SIMDIR,PARAMS,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',SIMDIR(4:end),PARAMS,'/']; BASIC.PARAMS = PARAMS; BASIC.SIMID = SIMID; BASIC.nrun = 1e8; BASIC.dt = DT; BASIC.tmax = TMAX; %time normalized to 1/omega_pe BASIC.maxruntime = str2num(CLUSTER.TIME(1:2))*3600 ... + str2num(CLUSTER.TIME(4:5))*60 ... + str2num(CLUSTER.TIME(7:8)); % Outputs parameters if RESTART; BASIC.RESTART = '.true.'; else; BASIC.RESTART = '.false.';end; OUTPUTS.nsave_0d = floor(1.0/SPS0D/DT); OUTPUTS.nsave_1d = -1; OUTPUTS.nsave_2d = floor(1.0/SPS2D/DT); OUTPUTS.nsave_3d = floor(1.0/SPS3D/DT); OUTPUTS.nsave_5d = floor(1.0/SPS5D/DT); OUTPUTS.nsave_cp = floor(1.0/SPSCP/DT); if W_DOUBLE; OUTPUTS.write_doubleprecision = '.true.'; else; OUTPUTS.write_doubleprecision = '.false.';end; if W_GAMMA; OUTPUTS.write_gamma = '.true.'; else; OUTPUTS.write_gamma = '.false.';end; if W_PHI; OUTPUTS.write_phi = '.true.'; else; OUTPUTS.write_phi = '.false.';end; if W_NA00; OUTPUTS.write_Na00 = '.true.'; else; OUTPUTS.write_Na00 = '.false.';end; if W_NAPJ; OUTPUTS.write_Napj = '.true.'; else; OUTPUTS.write_Napj = '.false.';end; if W_SAPJ; OUTPUTS.write_Sapj = '.true.'; else; OUTPUTS.write_Sapj = '.false.';end; if W_DENS; OUTPUTS.write_dens = '.true.'; else; OUTPUTS.write_dens = '.false.';end; if W_TEMP; OUTPUTS.write_temp = '.true.'; else; OUTPUTS.write_temp = '.false.';end; OUTPUTS.resfile0 = '''outputs'''; OUTPUTS.rstfile0 = '''checkpoint'''; OUTPUTS.job2load = JOB2LOAD; %% Create directories if ~exist(SIMDIR, 'dir') mkdir(SIMDIR) end if ~exist(BASIC.RESDIR, 'dir') mkdir(BASIC.RESDIR) end if ~exist(BASIC.MISCDIR, 'dir') mkdir(BASIC.MISCDIR) end %% Compile and WRITE input file INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC); nproc = 1; MAKE = 'cd ..; make; cd wk'; system(MAKE); %% disp(['Set up ',SIMID]); disp([resolution,gridname,degngrad]); if RESTART disp(['- restarting from JOBNUM = ',num2str(JOB2LOAD)]); else disp(['- starting from T = 0']); end diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m index bbaaf55..e77c9bf 100644 --- a/wk/analysis_3D.m +++ b/wk/analysis_3D.m @@ -1,566 +1,576 @@ addpath(genpath('../matlab')) % ... add for i_ = 1 % for ETA_ =[0.6:0.1:0.9] %% Load results if 1% Local results outfile =''; outfile =''; -outfile ='test_3D/100x50x4_L_60_P_2_J_1_eta_0.625_nu_1e-02_DGGK_CLOS_0_mu_2e-02'; -% outfile ='test_3D/100x50_L_60_P_2_J_1_eta_0.625_nu_1e-02_DGGK_CLOS_0_mu_2e-02'; +outfile =''; +outfile =''; +outfile =''; +outfile =''; +outfile =''; +outfile =''; +outfile =''; +outfile ='test_3D/100x50_L_100_P_2_J_1_eta_0.5_nu_1e-01_DGGK_mu_2e-01'; BASIC.RESDIR = ['../results/',outfile,'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/']; CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD); system(CMD); end if 0% Marconi results outfile =''; outfile =''; +outfile =''; +outfile =''; +outfile =''; +outfile =''; BASIC.RESDIR = ['../',outfile(46:end-8),'/']; BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/']; end %% JOBNUMMAX = 20; compile_results_3D %Compile the results from first output found to JOBNUMMAX if existing %% Retrieving max polynomial degree and sampling info Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe); Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi); Ns5D = numel(Ts5D); Ns3D = numel(Ts3D); % renaming and reshaping quantity of interest Ts5D = Ts5D'; Ts3D = Ts3D'; %% Build grids Nkx = numel(kx); Nky = numel(ky); [KY,KX] = meshgrid(ky,kx); Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky); dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1); KPERP2 = KY.^2+KX.^2; [~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky)); Lk = max(Lkx,Lky); dx = 2*pi/Lk; dy = 2*pi/Lk; Nx = max(Nkx,Nky); Ny = Nx; Nz = numel(z); x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x); y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y); [ZZ,RR] = meshgrid(y,x); %% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% disp('Analysis :') disp('- iFFT') % IFFT (Lower case = real space, upper case = frequency space) ne00 = zeros(Nx,Ny,Nz,Ns3D); % Gyrocenter density ni00 = zeros(Nx,Ny,Nz,Ns3D); dzTe = zeros(Nx,Ny,Nz,Ns3D); dzTi = zeros(Nx,Ny,Nz,Ns3D); dzni = zeros(Nx,Ny,Nz,Ns3D); np_i = zeros(Nx,Ny,Nz,Ns5D); % Ion particle density si00 = zeros(Nx,Ny,Nz,Ns5D); phi = zeros(Nx,Ny,Nz,Ns3D); dens_e = zeros(Nx,Ny,Nz,Ns3D); dens_i = zeros(Nx,Ny,Nz,Ns3D); temp_e = zeros(Nx,Ny,Nz,Ns3D); temp_i = zeros(Nx,Ny,Nz,Ns3D); drphi = zeros(Nx,Ny,Nz,Ns3D); dzphi = zeros(Nx,Ny,Nz,Ns3D); dr2phi = zeros(Nx,Ny,Nz,Ns3D); for it = 1:numel(Ts3D) for iz = 1:numel(z) NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); ne00(:,:,iz,it) = real(fftshift(ifft2((NE_),Nx,Ny))); ni00(:,:,iz,it) = real(fftshift(ifft2((NI_),Nx,Ny))); phi (:,:,iz,it) = real(fftshift(ifft2((PH_),Nx,Ny))); drphi(:,:,iz,it) = real(fftshift(ifft2(1i*KX.*(PH_),Nx,Ny))); dr2phi(:,:,iz,it)= real(fftshift(ifft2(-KX.^2.*(PH_),Nx,Ny))); dzphi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(PH_),Nx,Ny))); if(W_DENS && W_TEMP) DENS_E_ = DENS_E(:,:,iz,it); DENS_I_ = DENS_I(:,:,iz,it); TEMP_E_ = TEMP_E(:,:,iz,it); TEMP_I_ = TEMP_I(:,:,iz,it); dzni(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(DENS_I_),Nx,Ny))); dzTe(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_E_),Nx,Ny))); dzTi(:,:,iz,it) = real(fftshift(ifft2(1i*KY.*(TEMP_I_),Nx,Ny))); dens_e (:,:,iz,it) = real(fftshift(ifft2((DENS_E_),Nx,Ny))); dens_i (:,:,iz,it) = real(fftshift(ifft2((DENS_I_),Nx,Ny))); temp_e (:,:,iz,it) = real(fftshift(ifft2((TEMP_E_),Nx,Ny))); temp_i (:,:,iz,it) = real(fftshift(ifft2((TEMP_I_),Nx,Ny))); end end end Np_i = zeros(Nkx,Nky,Ns5D); % Ion particle density in Fourier space for it = 1:numel(Ts5D) [~, it2D] = min(abs(Ts3D-Ts5D(it))); Np_i(:,:,it) = 0; for ij = 1:Nji Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1)); Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it)); end np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nx,Ny))); end % Post processing disp('- post processing') % gyrocenter and particle flux from real space GFlux_xi = zeros(1,Ns3D); % Gyrocenter flux Gamma = GFlux_yi = zeros(1,Ns3D); % Gyrocenter flux Gamma = GFlux_xe = zeros(1,Ns3D); % Gyrocenter flux Gamma = GFlux_ye = zeros(1,Ns3D); % Gyrocenter flux Gamma = % Hermite energy spectrum epsilon_e_pj = zeros(Pe_max,Je_max,Ns5D); epsilon_i_pj = zeros(Pi_max,Ji_max,Ns5D); phi_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_maxx_avg = zeros(Nz,Ns3D); % Time evol. of the norm of phi phi_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi shear_maxx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_avgx_maxy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_maxx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear shear_avgx_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of shear Ne_norm = zeros(Pe_max,Je_max,Ns5D); % Time evol. of the norm of Napj Ni_norm = zeros(Pi_max,Ji_max,Ns5D); % . % Kperp spectrum interpolation %full kperp points kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]); % interpolated kperps nk_noAA = floor(2/3*numel(kx)); kp_ip = kx; [thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip); [xn,yn] = pol2cart(thg,rg); [ky_s, sortIdx] = sort(ky); [xc,yc] = meshgrid(ky_s,kx); phi_kp_t = zeros(numel(kp_ip),Nz,Ns3D); % for it = 1:numel(Ts3D) % Loop over 2D arrays for iz = 1:numel(z) NE_ = Ne00(:,:,iz,it); NI_ = Ni00(:,:,iz,it); PH_ = PHI(:,:,iz,it); phi_maxx_maxy(iz,it) = max( max(squeeze(phi(:,:,iz,it)))); phi_avgx_maxy(iz,it) = max(mean(squeeze(phi(:,:,iz,it)))); phi_maxx_avgy(iz,it) = mean( max(squeeze(phi(:,:,iz,it)))); phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it)))); shear_maxx_maxy(iz,it) = max( max(squeeze(-(dr2phi(:,:,iz,it))))); shear_avgx_maxy(iz,it) = max(mean(squeeze(-(dr2phi(:,:,iz,it))))); shear_maxx_avgy(iz,it) = mean( max(squeeze(-(dr2phi(:,:,iz,it))))); shear_avgx_avgy(iz,it) = mean(mean(squeeze(-(dr2phi(:,:,iz,it))))); GFlux_xi(iz,it) = sum(sum(ni00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly; GFlux_yi(iz,it) = sum(sum(-ni00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly; GFlux_xe(iz,it) = sum(sum(ne00(:,:,iz,it).*dzphi(:,:,iz,it)))*dx*dy/Lx/Ly; GFlux_ye(iz,it) = sum(sum(-ne00(:,:,iz,it).*drphi(:,:,iz,it)))*dx*dy/Lx/Ly; Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,iz,it))).^2,3)),xn,yn); phi_kp_t(:,iz,it) = mean(Z_rth,2); end end % for it = 1:numel(Ts5D) % Loop over 5D arrays [~, it2D] = min(abs(Ts3D-Ts5D(it))); Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkx/Nky; Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkx/Nky; epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4); epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4); end %% Compute primary instability growth rate disp('- growth rate') % Find max value of transport (end of linear mode) [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2); [~,itmax] = min(abs(Ts3D-tmax)); tstart = 0.1 * Ts3D(itmax); tend = 0.5 * Ts3D(itmax); [~,its3D_lin] = min(abs(Ts3D-tstart)); [~,ite3D_lin] = min(abs(Ts3D-tend)); g_I = zeros(Nkx,Nky,Nz); for ikx = 1:Nkx for iky = 1:Nky for iz = 1:Nz [g_I(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); end end end [gmax_I,ikmax_I] = max(max(g_I(1,:,:),[],2),[],3); kmax_I = abs(ky(ikmax_I)); Bohm_transport = ETAB/ETAN*gmax_I/kmax_I^2; %% Compute secondary instability growth rate disp('- growth rate') % Find max value of transport (end of linear mode) % [tmp,tmax] = max(GGAMMA_RI*(2*pi/Nx/Ny)^2); % [~,itmax] = min(abs(Ts2D-tmax)); % tstart = Ts2D(itmax); tend = 1.5*Ts2D(itmax); [~,its3D_lin] = min(abs(Ts3D-tstart)); [~,ite3D_lin] = min(abs(Ts3D-tend)); g_II = zeros(Nkx,Nky); for ikx = 1:Nkx for iky = 1 for iz = 1:Nz [g_II(ikx,iky,iz), ~] = LinearFit_s(Ts3D(its3D_lin:ite3D_lin),squeeze(abs(Ni00(ikx,iky,iz,its3D_lin:ite3D_lin)))); end end end [gmax_II,ikmax_II] = max(max(g_II(1,:,:),[],2),[],3); kmax_II = abs(kx(ikmax_II)); %% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default_plots_options disp('Plots') FMT = '.fig'; if 1 %% Time evolutions and growth rate fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS]; set(gcf, 'Position', [100, 100, 900, 800]) subplot(111); suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... ' $\mu_{hd}=$',num2str(MU)]); subplot(421); for ip = 1:Pe_max for ij = 1:Je_max plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_e^{',num2str(ip-1),num2str(ij-1),'}$']; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,... 'Color',clr,'LineStyle',lstyle{1}); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$'); subplot(423) for ip = 1:Pi_max for ij = 1:Ji_max plt = @(x) squeeze(x(ip,ij,:)); plotname = ['$N_i^{',num2str(ip-1),num2str(ij-1),'}$']; clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,... 'Color',clr,'LineStyle',lstyle{1}); hold on; end end grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$') subplot(222) plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on; plot(Ts0D,PGAMMA_RI*(2*pi/Nx/Ny)^2); legend(['Gyro. flux';'Part. flux']); grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$') if(1) subplot(223) plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on; plot(kx,max(g_II(:,1,:),[],3),'x-','DisplayName','Second. instability'); hold on; plot([max(ky)*2/3,max(ky)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA'); grid on; xlabel('$k\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show'); ylim([0,max(max(g_I(1,:,:)))]); xlim([0,max(ky)]); shearplot = 426; phiplot = 428; else shearplot = 223; phiplot = 224; end subplot(shearplot) plt = @(x) mean(x,1); clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s)$'); hold on; plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s \rangle_z$'); hold on; plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s \rangle_r$'); hold on; plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s \rangle_{r,z}$'); hold on; grid on; xlabel('$t c_s/R$'); ylabel('$shear$'); subplot(phiplot) clr = line_colors(min(ip,numel(line_colors(:,1))),:); lstyle = line_styles(min(ij,numel(line_styles))); plot(Ts3D,plt(phi_maxx_maxy),'DisplayName','$\max_{r,z}(\phi)$'); hold on; plot(Ts3D,plt(phi_maxx_avg),'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on; plot(Ts3D,plt(phi_avgx_maxy),'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on; plot(Ts3D,plt(phi_avgx_avgy),'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on; grid on; xlabel('$t c_s/R$'); ylabel('$E.S. pot$'); save_figure end if 1 %% Space time diagramm (fig 11 Ivanov 2020) TAVG = 5000; % Averaging time duration %Compute steady radial transport tend = Ts0D(end); tstart = tend - TAVG; [~,its0D] = min(abs(Ts0D-tstart)); [~,ite0D] = min(abs(Ts0D-tend)); SCALE = (2*pi/Nx/Ny)^2; gamma_infty_avg = mean(PGAMMA_RI(its0D:ite0D))*SCALE; gamma_infty_std = std (PGAMMA_RI(its0D:ite0D))*SCALE; % Compute steady shearing rate tend = Ts3D(end); tstart = tend - TAVG; [~,its2D] = min(abs(Ts3D-tstart)); [~,ite2D] = min(abs(Ts3D-tend)); shear_infty_avg = mean(mean(shear_maxx_avgy(:,its2D:ite2D),1)); shear_infty_std = std (mean(shear_maxx_avgy(:,its2D:ite2D),1)); % plots fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) subplot(311) % yyaxis left plot(Ts0D,PGAMMA_RI*SCALE,'DisplayName','$\langle n_i d\phi/dz \rangle_z$'); hold on; plot(Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',... 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]); grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_r$') ylim([0,5*abs(gamma_infty_avg)]); xlim([0,Ts0D(end)]); title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... ' $\mu_{hd}=$',num2str(MU)]); % subplot(312) clr = line_colors(1,:); lstyle = line_styles(1); % plt = @(x_) mean(x_,1); plt = @(x_) x_(1,:); plot(Ts3D,plt(shear_maxx_maxy),'DisplayName','$\max_{r,z}(s_\phi)$'); hold on; plot(Ts3D,plt(shear_maxx_avgy),'DisplayName','$\max_{r}\langle s_\phi\rangle_z$'); hold on; plot(Ts3D,plt(shear_avgx_maxy),'DisplayName','$\max_{z}\langle s_\phi\rangle_r$'); hold on; plot(Ts3D,plt(shear_avgx_avgy),'DisplayName','$\langle s_\phi\rangle_{r,z}$'); hold on; plot(Ts3D(its2D:ite2D),ones(ite2D-its2D+1,1)*shear_infty_avg, '-k',... 'DisplayName',['$s^{\infty} = $',num2str(shear_infty_avg),'$\pm$',num2str(shear_infty_std)]); ylim([0,shear_infty_avg*5.0]); xlim([0,Ts0D(end)]); grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show'); subplot(313) [TY,TX] = meshgrid(x,Ts3D); pclr = pcolor(TX,TY,squeeze((mean(dr2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_r^2\phi\rangle_z$)') %colorbar; caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$r/\rho_s$'); colormap(bluewhitered(256)) save_figure end if 1 %% Space time diagramms tstart = 0; tend = Ts3D(end); [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); trange = itstart:itend; [TY,TX] = meshgrid(x,Ts3D(trange)); fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600]) subplot(211) % pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; shading interp colormap hot; caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]); caxis([0.0,0.01]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z'; xticks([]); ylabel('$r/\rho_s$') % legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$') title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... ', $L=',num2str(L),'$, $N=',num2str(Nx),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... ' $\mu_{hd}=$',num2str(MU)]); subplot(212) pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,1,trange),2))'); set(pclr, 'edgecolor','none'); colorbar; fieldmax = max(max(mean(abs(drphi(:,:,1,its2D:ite2D)),2))); caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z'; xlabel('$t c_s/R$'), ylabel('$r/\rho_s$') % legend('Zonal flow $\langle \partial_r\phi\rangle_z$') save_figure end if 0 %% Averaged shear and Reynold stress profiles trange = its2D:ite2D; % trange = 100:200; figure; plt = @(x) squeeze(mean(mean(real(x(:,:,trange)),2),3))/max(abs(squeeze(mean(mean(real(x(:,:,trange)),2),3)))); plot(r,plt(-dr2phi),'-k','DisplayName','Zonal shear'); hold on; plot(r,plt(-drphi.*dzphi),'-','DisplayName','$\Pi_\phi$'); hold on; % plot(r,plt(-drphi.*dzTe),'-','DisplayName','$\Pi_{Te}$'); hold on; plot(r,plt(-drphi.*dzTi),'-','DisplayName','$\Pi_{Ti}$'); hold on; plot(r,plt(-drphi.*dzphi-drphi.*dzTi),'-','DisplayName','$\Pi_\phi+\Pi_{Ti}$'); hold on; % plot(r,plt(-drphi.*dzphi-drphi.*dzTi-drphi.*dzTe),'-','DisplayName','$\Pi_\phi+\Pi_{Te}+\Pi_{Ti}$'); hold on; xlim([-L/2,L/2]); xlabel('$r/\rho_s$'); grid on; legend('show') end if 1 %% |phi_k|^2 spectra (Kobayashi 2015 fig 3) tstart = 0.8*Ts3D(end); tend = Ts3D(end); [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); trange = itstart:itend; %full kperp points phi_k_2 = reshape(mean(mean(abs(PHI(:,:,:,trange)),3),4).^2,[numel(KX),1]); kperp = reshape(sqrt(KX.^2+KY.^2),[numel(KX),1]); % interpolated kperps nk_noAA = floor(2/3*numel(kx)); kp_ip = kx; [thg, rg] = meshgrid(linspace(0,pi,2*nk_noAA),kp_ip); [xn,yn] = pol2cart(thg,rg); [ky_s, sortIdx] = sort(ky); [xc,yc] = meshgrid(ky_s,kx); Z_rth = interp2(xc,yc,squeeze(mean((abs(PHI(:,sortIdx,trange))).^2,3)),xn,yn); phi_kp = mean(Z_rth,2); Z_rth = interp2(xc,yc,squeeze(mean((abs(Ni00(:,sortIdx,trange))).^2,3)),xn,yn); ni00_kp = mean(Z_rth,2); Z_rth = interp2(xc,yc,squeeze(mean((abs(Ne00(:,sortIdx,trange))).^2,3)),xn,yn); ne00_kp = mean(Z_rth,2); %for theorical trends a1 = phi_kp(2)*kp_ip(2).^(13/3); a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2); fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500]) % scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on; plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on; plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on; plot(kp_ip,ne00_kp, '^','DisplayName','$\langle|N_{e,k}^{00}|^2\rangle_{k_\perp}$'); hold on; plot(kp_ip,a1*kp_ip.^(-13/3),'-','DisplayName','$k^{-13/3}$'); plot(kp_ip,a2/100*kp_ip.^(-3)./(1+kp_ip.^2).^2,'-','DisplayName','$k^{-3}/(1+k^2)^2$'); set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest') title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta=$',num2str(ETAB/ETAN),... ', $L=',num2str(L),'$, $N=',num2str(Nx),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',... ' $\mu_{hd}=$',num2str(MU)]}); xlim([0.1,10]); save_figure clear Z_rth phi_kp ni_kp Ti_kp end %% t0 =00; [~, it02D] = min(abs(Ts3D-t0)); [~, it05D] = min(abs(Ts5D-t0)); skip_ = 1; DELAY = 1e-3*skip_; -FRAMES_2D = it02D:skip_:numel(Ts3D); +FRAMES_3D = it02D:skip_:numel(Ts3D); FRAMES_5D = it05D:skip_:numel(Ts5D); %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0 %% part density electron GIFNAME = ['ne',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(mean(dens_e,3)); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; -% FIELD = real(dens_e(:,:,2,:)); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(mean(dens_e,3)); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; +% FIELD = real(dens_e(:,:,2,:)); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_3D; FIELDNAME = '$n_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% part temperature electron GIFNAME = ['Te',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(temp_e); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(temp_e); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$T_e$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% part density ion GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(dens_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(dens_i); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% part temperature ion GIFNAME = ['Ti',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(temp_i); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(temp_i); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$T_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% GC Density ion GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(ni00); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$n_i^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% GC Density electrons GIFNAME = ['ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; -FIELD = real(ne00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(ne00); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$n_e^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; % create_gif create_mov end if 0 %% Phi real space GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0; -FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = real(phi); X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif % create_mov end if 0 %% shear GIFNAME = ['shear_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1; -FIELD = -dr2phi; X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D; +FIELD = -dr2phi; X = RR; Y = ZZ; T = Ts3D; FRAMES = FRAMES_3D; FIELDNAME = '$s$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$'; create_gif end if 0 %% phi averaged on z GIFNAME = ['phi_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1; -FIELD =(squeeze(mean(real(phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D; -X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); +FIELD =(squeeze(mean(real(phi),2))); linestyle = '-.k'; FRAMES = FRAMES_3D; +X = (r); T = Ts3D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$'; create_gif_1D_phi end if 0 %% flow averaged on z GIFNAME = ['zf_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1; -FIELD =(squeeze(mean(real(drphi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D; -X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); +FIELD =(squeeze(mean(real(drphi),2))); linestyle = '-.k'; FRAMES = FRAMES_3D; +X = (r); T = Ts3D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$'; create_gif_1D_phi end if 0 %% shear averaged on z GIFNAME = ['shear_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1; -FIELD =(squeeze(mean(real(dr2phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D; -X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); +FIELD =(squeeze(mean(real(dr2phi),2))); linestyle = '-.k'; FRAMES = FRAMES_3D; +X = (r); T = Ts3D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r); FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$'; create_gif_1D_phi end if 0 %% phi kperp spectrum GIFNAME = ['phi2_kp',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0; -FIELD =log10(phi_kp_t); linestyle = '-'; FRAMES = FRAMES_2D; -X = kp_ip; T = Ts2D; YMIN = -20; YMAX = 10; XMIN = min(kx); XMAX = max(kx); +FIELD =log10(phi_kp_t); linestyle = '-'; FRAMES = FRAMES_3D; +X = kp_ip; T = Ts3D; YMIN = -20; YMAX = 10; XMIN = min(kx); XMAX = max(kx); FIELDNAME = '$\log_{10}|\tilde\phi_k|^2$'; XNAME = '$k_r\rho_s$'; create_gif_1D end if 0 %% Density ion frequency -GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D; -FIELD =ifftshift((abs(Ni00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts2D; +GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_3D; +FIELD =ifftshift((abs(Ni00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts3D; FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% Density electron frequency -GIFNAME = ['Ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D; -FIELD =ifftshift((abs(Ne00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts2D; +GIFNAME = ['Ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_3D; +FIELD =ifftshift((abs(Ne00)),2); X = fftshift(kx,2); Y = fftshift(ky,2); T = Ts3D; FIELDNAME = '$N_e^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$'; create_gif end if 0 %% kx vs P Si GIFNAME = ['Sip0_kx',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) squeeze(max((abs(x)),[],4)); FIELD =plt(Sipj(:,1,:,:,:)); X = kx'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$, $\sum_z$'; create_gif_imagesc end if 0 %% maxky, kx vs p, for all Nipj over time GIFNAME = ['Nipj_kx',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) squeeze(sum((abs(x)),4)); FIELD = plt(Nipj); X = kx'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$'; create_gif_5D end if 0 %% maxkx, ky vs p, for all Nipj over time GIFNAME = ['Nipj_ky',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; plt = @(x) fftshift(squeeze(sum((abs(x)),3)),3); FIELD = plt(Nipj); X = sort(ky'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D; FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, $\sum_r$'; create_gif_5D end %% % ZF_fourier_analysis end \ No newline at end of file diff --git a/wk/local_run.m b/wk/local_run.m index 76e0198..d03b5f8 100644 --- a/wk/local_run.m +++ b/wk/local_run.m @@ -1,74 +1,74 @@ addpath(genpath('../matlab')) % ... add %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS -NU = 1e-2; % Collision frequency -ETAN = 1.6; % Density gradient drive (R/Ln) +NU = 1e-1; % Collision frequency +ETAN = 2.0; % Density gradient drive (R/Ln) NU_HYP = 1.0; %% GRID PARAMETERS N = 100; % Frequency gridpoints (Nkx = N/2) -L = 60; % Size of the squared frequency domain -Nz = 1; % number of perpendicular planes (parallel grid) -Lz = 1; % size of the parallel length +L = 100; % Size of the squared frequency domain +Nz = 10; % number of perpendicular planes (parallel grid) +Lz = 2*pi; % size of the parallel length P = 2; J = 1; MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f %% TIME PARAMETERS TMAX = 100; % Maximal time unit -DT = 2e-2; % Time step +DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for profiler SPS2D = 1; % Sampling per time unit for 2D arrays SPS3D = 1; % Sampling per time unit for 3D arrays SPS5D = 1/200; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints/10 RESTART = 0; % To restart from last checkpoint JOB2LOAD= 0; %% OPTIONS AND NAMING % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK) CO = 1; CLOS = 0; % Closure model (0: =0 truncation) NL_CLOS = 0; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS) SIMID = 'test_3D'; % Name of the simulation % SIMID = 'v2.8_barnes'; % Name of the simulation % SIMID = ['v2.8_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1) INIT_ZF = 0; ZF_AMP = 0.0; %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 1; W_PHI = 1; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; W_DENS = 1; W_TEMP = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% unused PMAXE = P; % Highest electron Hermite polynomial degree JMAXE = J; % Highest '' Laguerre '' PMAXI = P; % Highest ion Hermite polynomial degree JMAXI = J; % Highest '' Laguerre '' KERN = 0; % Kernel model (0 : GK) KX0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst. KXEQ0 = 0; % put kx = 0 KPAR = 0.0; % Parellel wave vector component LAMBDAD = 0.0; kmax = N*pi/L;% Highest fourier mode HD_CO = 0.5; % Hyper diffusivity cutoff ratio % kmaxcut = 2.5; MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient NOISE0 = 1.0e-5; TAU = 1.0; % e/i temperature ratio ETAT = 0.0; % Temperature gradient ETAB = 1.0; % Magnetic gradient (1.0 to set R=LB) INIT_PHI= 1; % Start simulation with a noisy phi and moments %% Setup and file management setup system('rm fort.90'); outfile = [BASIC.RESDIR,'out.txt']; disp(outfile);