diff --git a/matlab/setup.m b/matlab/setup.m index 821068b..07937e1 100644 --- a/matlab/setup.m +++ b/matlab/setup.m @@ -1,149 +1,150 @@ %% ________________________________________________________________________ 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 resolution GRID.q0 = q0; % q factor GRID.shear = shear; % shear GRID.eps = eps; % inverse aspect ratio % 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; if (WIPE_TURB); INITIAL.wipe_turb = '.true.'; else; INITIAL.wipe_turb = '.false.';end; if (INIT_BLOB); INITIAL.init_blob = '.true.'; else; INITIAL.init_blob = '.false.';end; 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 - INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5''']; + INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5''']; +% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5''']; elseif (CO == -1) % DGDK disp('Warning, DGDK not debugged') end % Naming and creating input file switch abs(CO) case 0; CONAME = 'LB'; case 1; CONAME = 'DG'; case 2; CONAME = 'SG'; case 3; CONAME = 'PA'; case 4; CONAME = 'FC'; otherwise; CONAME ='UK'; end if (CO <= 0); CONAME = [CONAME,'DK']; else; CONAME = [CONAME,'GK']; 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,'_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),'_q0_',num2str(q0),'_']; end 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/linear_study.m b/wk/linear_study.m index 81d5000..c0ac260 100644 --- a/wk/linear_study.m +++ b/wk/linear_study.m @@ -1,215 +1,224 @@ -for NU = [0.1] -for ETAN = [2.0] -for CO = [1] %clear all; addpath(genpath('../matlab')) % ... add default_plots_options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Set Up parameters CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PHYSICAL PARAMETERS % NU = 1.0; % Collision frequency TAU = 1.0; % e/i temperature ratio ETAB = 1.0; % ETAN = ETAN; % Density gradient ETAT = 0.0; % Temperature gradient NU_HYP = 0.0; % Hyperdiffusivity coefficient LAMBDAD = 0.0; NOISE0 = 1.0e-5; %% GRID PARAMETERS N = 50; % Frequency gridpoints (Nkx = N/2) -L = 100; % Size of the squared frequency domain +L = 75; % Size of the squared frequency domain KXEQ0 = 1; % put kx = 0 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 PARMETERS TMAX = 200; % Maximal time unit DT = 1e-2; % Time step SPS0D = 1; % Sampling per time unit for 2D arrays SPS2D = 1; % Sampling per time unit for 2D arrays SPS5D = 1/50; % Sampling per time unit for 5D arrays SPSCP = 0; % Sampling per time unit for checkpoints RESTART = 0; % To restart from last checkpoint JOB2LOAD= 00; %% OPTIONS SIMID = 'v3.1_lin_analysis'; % Name of the simulation NON_LIN = 0 *(1-KXEQ0); % activate non-linearity (is cancelled if KXEQ0 = 1) % Collision operator % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK) % CO = 2; INIT_ZF = 0; ZF_AMP = 0.0; CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P) NL_CLOS = 0; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS) KERN = 0; % Kernel model (0 : GK) INIT_PHI= 0; % Start simulation with a noisy phi INIT_ZF = 0; % Start simulation with a noisy phi %% OUTPUTS W_DOUBLE = 0; W_GAMMA = 0; W_PHI = 0; W_NA00 = 1; W_NAPJ = 1; W_SAPJ = 0; W_DENS = 0; W_TEMP = 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % unused % DK = 0; % Drift kinetic model (put every kernel_n to 0 except n=0 to 1) JOBNUM = 00; KPAR = 0.0; % Parellel wave vector component HD_CO = 0.5; % Hyper diffusivity cutoff ratio kmax = N*pi/L;% Highest fourier mode MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient Nz = 1; % number of perpendicular planes (parallel grid) q0 = 1.0; % safety factor shear = 0.0; % magnetic shear eps = 0.0; % inverse aspect ratio %% PARAMETER SCANS -if 1 + +if 0 %% Parameter scan over PJ % PA = [2 4 6 10]; % JA = [1 2 3 5]; PA = [6]; JA = [3]; DTA= DT*ones(size(JA));%./sqrt(JA); % DTA= DT; mup_ = MU_P; muj_ = MU_J; Nparam = numel(PA); param_name = 'PJ'; gamma_Ni00 = zeros(Nparam,floor(N/2)+1); gamma_Nipj = zeros(Nparam,floor(N/2)+1); Bohm_transport = zeros(Nparam,1); Ni00_ST = zeros(Nparam,floor(N/2)+1,floor(SPS2D*TMAX)); for i = 1:Nparam % Change scan parameter PMAXE = PA(i); PMAXI = PA(i); JMAXE = JA(i); JMAXI = JA(i); DT = DTA(i); setup % Run linear simulation % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz_3 1 1; cd ../../../wk']) system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']) % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 2 3; cd ../../../wk']) % Load and process results %% filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results tend = Ts3D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts3D-tstart)); [~,itend] = min(abs(Ts3D-tend)); for ikx = 1:N/2+1 gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend)))))); Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:))); end tend = Ts5D(end); tstart = 0.4*tend; [~,itstart] = min(abs(Ts5D-tstart)); [~,itend] = min(abs(Ts5D-tend)); for ikx = 1:N/2+1 gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2))); end gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0)); gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0)); % kymax = abs(kx(ikymax)); % Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; % Clean output system(['rm -r ',BASIC.RESDIR]); end if 1 %% Plot SCALE = 1;%sqrt(2); fig = figure; FIGNAME = 'linear_study'; plt = @(x) x; % subplot(211) for i = 1:Nparam clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),... 'Color',clr,... 'LineStyle',linestyle{1},'Marker','^',... 'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']); hold on; end grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]); title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$']) legend('show'); xlim([0.01,10]) saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']); saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']); end end -if 0 -%% Parameter scan over eta_B +if 1 +%% Parameter scan over CO PMAXE = 6; PMAXI = 6; JMAXE = 3; JMAXI = 3; -TMAX = 400; -DT = 0.001; -eta_B = [0.45, 0.5, 0.55, 0.6, 0.65]; -Nparam = numel(eta_B); -param_name = 'etaB'; -CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) -NU = 1e-2; % Collision frequency +TMAX = 200; +DT = 0.01; +CO_A = [4]; +CONAME_A = {}; +Nparam = numel(CO_A); +param_name = 'CO'; +ETAN = 2.0; +NU = 1e-1; % Collision frequency Bohm_transport = zeros(Nparam,1); -gamma_Ni = zeros(Nparam,N); +gamma_Ni00 = zeros(Nparam,floor(N/2)+1); +gamma_Nipj = zeros(Nparam,floor(N/2)+1); for i = 1:Nparam % Change scan parameter - ETAB = eta_B(i); + CO = CO_A(i); setup + CONAME_A{i} = CONAME; % Run linear simulation system(... - ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ./../../../bin/helaz; cd ../../../wk']... + ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk']... ) - % Load and process results + %% Load an process results + filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5']; load_results - tend = Ts2D(end); tstart = 0.4*tend; + tend = Ts3D(end); tstart = 0.4*tend; + [~,itstart] = min(abs(Ts3D-tstart)); + [~,itend] = min(abs(Ts3D-tend)); for ikx = 1:N/2+1 - gamma_Ni(i,ikx) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikx,1,:))),tstart,tend); - Ni00_ST(i,ikx,1:numel(Ts2D)) = squeeze((Ni00(ikx,1,:))); + gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend)))))); + Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:))); end - gamma_Ni(i,:) = real(gamma_Ni(i,:) .* (gamma_Ni(i,:)>=0.0)); - [gmax,ikymax] = max(gamma_Ni(i,:)); - kymax = abs(kx(ikymax)); - Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; + tend = Ts5D(end); tstart = 0.4*tend; + [~,itstart] = min(abs(Ts5D-tstart)); + [~,itend] = min(abs(Ts5D-tend)); + for ikx = 1:N/2+1 + gamma_Nipj(i,ikx) = LinearFit_s(Ts5D(itstart:itend)',squeeze(max(max(abs(Nipj(:,:,ikx,1,itstart:itend)),[],1),[],2))); + end + gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0)); + gamma_Nipj(i,:) = real(gamma_Nipj(i,:) .* (gamma_Nipj(i,:)>=0.0)); +% kymax = abs(kx(ikymax)); +% Bohm_transport(i) = ETAB/ETAN*gmax/kymax^2; % Clean output - system(['rm -r ',BASIC.RESDIR]) + system(['rm -r ',BASIC.RESDIR]); end -if 0 +if 1 %% Plot +SCALE = 1;%sqrt(2); fig = figure; FIGNAME = 'linear_study'; -plt = @(x) circshift(x,N/2-1); -for i = 1:Nparam - clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); - linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); - plot(plt(ky),plt(gamma_Ni(i,:)),... - 'Color',clr,... - 'LineStyle',linestyle{1},... - 'DisplayName',['$\eta_B=$',num2str(eta_B(i))]); - hold on; -end -grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma(N_i^{00})\rho_2/c_s$'); xlim([0.0,max(ky)]); -title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... - ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) -legend('show') +plt = @(x) x; +% subplot(211) + for i = 1:Nparam + clr = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:); + linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1); + semilogx(plt(SCALE*kx(2:numel(kx))),plt(gamma_Ni00(i,2:end)),... + 'Color',clr,... + 'LineStyle',linestyle{1},'Marker','^',... + 'DisplayName',[CONAME_A{i}]); + hold on; + end + grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})L_\perp/c_s$'); xlim([0.0,max(kx)]); + title(['$\eta=',num2str(ETAB/ETAN),'$, $\nu=',num2str(NU),'$',', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']) + legend('show'); xlim([0.01,10]) saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']); +saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']); end -if 1 +if 0 %% Plot fig = figure; FIGNAME = 'mixing_length'; plot(eta_B, Bohm_transport) grid on; xlabel('$L_n/L_B$'); ylabel('$\eta\gamma_{max}/k_{max}^2$'); title(['$P_e=',num2str(PMAXE),'$',', $J_e=',num2str(JMAXE),'$',... ', $P_i=',num2str(PMAXE),'$',', $J_i=',num2str(JMAXI),'$']) saveas(fig,[SIMDIR,FIGNAME,'_vs_',param_name,'_',PARAMS,'.fig']); end %% -end -end -end end \ No newline at end of file diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m index a6f583a..b81ce09 100644 --- a/wk/plot_cosol_mat.m +++ b/wk/plot_cosol_mat.m @@ -1,89 +1,126 @@ % MAT = results.iCa; % figure % suptitle('FCGK,P=6,J=3'); % subplot(221) % imagesc(log(abs(MAT))); % title('log abs') % subplot(222) % imagesc(imag(MAT)); colormap(bluewhitered) % title('imag'); colorbar % subplot(223) % imagesc(imag(MAT)>0); % title('imag$>$0'); % subplot(224) % imagesc(imag(MAT)<0); % title('imag$<$0'); % % %% SGGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'; SGGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj'); SGGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); SGGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=0'); % MAT = 1i*SGGK_ei; suptitle('SGGK ei,P=20,J=10'); % MAT = 1i*SGGK_ie; suptitle('SGGK ie,P=20,J=10'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); %% PAGK P_ = 20; J_ = 10; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'; PAGK_self = h5read(mat_file_name,'/00000/Caapj/Ceepj'); % PAGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); % PAGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=0'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); - %% SGGK + %% FCGK P_ = 6; J_ = 3; mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5'; -FCGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj'); -FCGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); -FCGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); - +kp = 1.2; matidx = round(kp/(8/150)); +matidx = sprintf('%5.5i',matidx); +FCGK_self = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); +FCGK_ei = h5read(mat_file_name,['/',matidx,'/Ceipj/CeipjF'])+h5read(mat_file_name,'/00000/Ceipj/CeipjT'); +FCGK_ie = h5read(mat_file_name,['/',matidx,'/Ciepj/CiepjF'])+h5read(mat_file_name,'/00000/Ciepj/CiepjT'); figure -MAT = 1i*FCGK_self; suptitle('FCGK ii,P=20,J=10, k=0'); -% MAT = 1i*SGGK_ei; suptitle('FCGK ei,P=20,J=10'); -% MAT = 1i*SGGK_ie; suptitle('FCGK ie,P=20,J=10'); +MAT = 1i*FCGK_self; suptitle(['FCGK ii,P=20,J=10, k=',num2str(kp)]); +% MAT = 1i*FCGK_ei; suptitle('FCGK ei,P=20,J=10'); +% MAT = 1i*FCGK_ie; suptitle('FCGK ie,P=20,J=10'); subplot(221) imagesc(abs(MAT)); title('log abs') subplot(222) imagesc(imag(MAT)); colormap(bluewhitered) title('imag'); colorbar subplot(223) imagesc(imag(MAT)>0); title('imag$>$0'); subplot(224) imagesc(imag(MAT)<0); title('imag$<$0'); - \ No newline at end of file + +%% Eigenvalue analysis +if 0 +%% +mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',... + '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',... + '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_4.h5',... + '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5',... + }; +CONAME_A = {'SG','PA','FC NFLR 4', 'FC NFLR k2'}; +kp_a = linspace(0,8,149); +gmax = zeros(size(kp_a)); +wmax = zeros(size(kp_a)); +figure +for j_ = 1:4 + mat_file_name = mfns{j_}; + i = 1; + for kp = kp_a + matidx = floor(kp/(8/149)); + matidx = sprintf('%5.5i',matidx);disp(matidx) + MAT = h5read(mat_file_name,['/',matidx,'/Caapj/Ciipj']); + gmax(i) = max(real(eig(MAT))); + wmax(i) = max(imag(eig(MAT))); + i = i + 1; + end + subplot(121) + plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on; + subplot(122) + plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on; +end + subplot(121) +legend('show'); grid on; +xlabel('$k_\perp$'); ylabel('$\gamma_{max}$ from Eig(iCa)') + subplot(122) +legend('show'); grid on; +xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)') +end