diff --git a/Makefile b/Makefile
index 54fd33f091d4b391f8559a102801e31820397a33..8d514dcc7d7318e0c35b61d39927c75427eb40af 100644
--- a/Makefile
+++ b/Makefile
@@ -60,7 +60,7 @@ $(OBJDIR)/utility_mod.o
$(EXEC): $(FOBJ)
$(F90) $(LDFLAGS) $(OBJDIR)/*.o $(EXTMOD) $(EXTINC) $(EXTLIBS) -o $@
- $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
+ $(OBJDIR)/advance_field.o : src/advance_field.F90 $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/advance_field.F90 -o $@
$(OBJDIR)/array_mod.o : src/array_mod.F90 $(OBJDIR)/prec_const_mod.o
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 299d65cd32c4aa99f6d01735ea0f8f34f2c07650..0f49d65cc61f09e2d28959743164b02ddad1bcfb 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -7,9 +7,9 @@ YNAME = latexize(YNAME);
FIELDNAME = latexize(FIELDNAME);
% Set colormap boundaries
-hmax = max(max(max(FIELD)));
-hmin = min(min(min(FIELD)));
-
+hmax = max(max(max(FIELD(:,:,FRAMES))));
+hmin = min(min(min(FIELD(:,:,FRAMES))));
+scale = -1;
flag = 0;
if hmax == hmin
disp('Warning : h = hmin = hmax = const')
@@ -25,17 +25,28 @@ fig = figure('Color','white','Position', [100, 100, 400, 400]);
in = 1;
nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
for n = FRAMES % loop over selected frames
- scale = max(max(abs(FIELD(:,:,n))));
- pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot
+ scale = max(max(abs(FIELD(:,:,n)))); % Scaling to normalize
+ if NORMALIZED
+ pclr = pcolor(X,Y,FIELD(:,:,n)/scale); % frame plot\
+ caxis([-1,1]);
+ else
+ pclr = pcolor(X,Y,FIELD(:,:,n)); % frame plot\
+ if CONST_CMAP
+ caxis([-1,1]*max(abs([hmin hmax]))); % adaptive color map
+ else
+ caxis([-1,1]*scale); % adaptive color map
+ end
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scale = ',sprintf('%.1e',scale)]);
+ end
+ title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
+ ,', scaling = ',sprintf('%.1e',scale)]);
if INTERP
shading interp;
end
set(pclr, 'edgecolor','none'); axis square;
- caxis([-1,1]);
colormap(bluewhitered)
xlabel(XNAME); ylabel(YNAME); %colorbar;
- title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale)]);
drawnow
% Capture the plot as an image
frame = getframe(fig);
diff --git a/matlab/load_params.m b/matlab/load_params.m
index 1281941854829eb4bfb2be00835923c667b0951f..1e93d9e76e25a517dda6ea5777f661c7a48ef3cc 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -27,13 +27,16 @@ else
NON_LIN = 0;
end
-if (CO == -3); CONAME = 'PADK';
-elseif(CO == -2); CONAME = 'SGDK';
-elseif(CO == -1); CONAME = 'DGDK';
-elseif(CO == 0); CONAME = 'LBGK';
-elseif(CO == 1); CONAME = 'DGGK';
-elseif(CO == 2); CONAME = 'SGGK';
-elseif(CO == 3); CONAME = 'PAGK';
+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.';
diff --git a/matlab/plots/plot_radial_transport_and_shear.m b/matlab/plots/plot_radial_transport_and_shear.m
index 30a3128e1ad8cc04b1a850a75bbf9ef45e380e25..86bd431911cf1e85eaf6bf1c80930ea015efb5de 100644
--- a/matlab/plots/plot_radial_transport_and_shear.m
+++ b/matlab/plots/plot_radial_transport_and_shear.m
@@ -39,6 +39,7 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
subplot(313)
[TY,TX] = meshgrid(x,Ts3D);
- pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('Shear ($\langle \partial_x^2\phi\rangle_y$)') %colorbar;
- caxis(1*shear_infty_avg*[-1 1]); xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
+ pclr = pcolor(TX,TY,squeeze((mean(dx2phi(:,:,1,:),2)))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_x^2\phi\rangle_y$') %colorbar;
+ caxis([-1 1]);
+ xlabel('$t c_s/R$'), ylabel('$x/\rho_s$'); colormap(bluewhitered(256))
save_figure
\ No newline at end of file
diff --git a/matlab/plots/plot_space_time_diagrams.m b/matlab/plots/plot_space_time_diagrams.m
index 7f101ff0d2809e9b12bf4704981a679e0e9b6868..4fb784998788f35fb1aa13ef5666851f51cd9477 100644
--- a/matlab/plots/plot_space_time_diagrams.m
+++ b/matlab/plots/plot_space_time_diagrams.m
@@ -21,6 +21,7 @@ fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 10
'DisplayName','$\langle \partial_r\phi\rangle_z$'); colorbar;
fieldmax = max(max(mean(abs(dxphi(:,:,1,its2D:ite2D)),2)));
caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle v_{E\times B,z}\rangle_z';
+ colormap(bluewhitered)
xlabel('$t c_s/R$'), ylabel('$x/\rho_s$')
% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
save_figure
\ No newline at end of file
diff --git a/matlab/plots/plot_time_evolution_and_gr.m b/matlab/plots/plot_time_evolution_and_gr.m
index e016520b4e1807191c8a844fed3bbc30832531a7..4427e3f250c2d32bfeb0d892169a9a9ccf1b98ad 100644
--- a/matlab/plots/plot_time_evolution_and_gr.m
+++ b/matlab/plots/plot_time_evolution_and_gr.m
@@ -29,10 +29,12 @@ subplot(111);
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}$')
+% plot(Ts0D,GGAMMA_RI*(2*pi/Nx/Ny)^2); hold on;
+ yyaxis left; ylabel('$\Gamma_x$');
+ plot(Ts3D,squeeze(sum(sum(sum(Gamma_x,1),2),3)));
+ yyaxis right; ylabel('$Q_x$');
+ plot(Ts3D,squeeze(sum(sum(sum(Q_x,1),2),3)));
+ grid on; xlabel('$t c_s/R$');
if(~isnan(max(max(g_I(1,:,:)))))
subplot(223)
plot(ky,max(g_I(1,:,:),[],3),'-','DisplayName','Primar. instability'); hold on;
diff --git a/matlab/post_processing.m b/matlab/post_processing.m
index 5b50bd397eca0044f1e79c309dde084f008d39ac..1a74cb29e04777d1391abd142da91d5887f90484 100644
--- a/matlab/post_processing.m
+++ b/matlab/post_processing.m
@@ -86,7 +86,13 @@ end
disp('- post processing')
% particle flux
Gamma_x= zeros(Nx,Ny,Nz,Ns3D); % Radial particle transport
+Gamma_y= zeros(Nx,Ny,Nz,Ns3D); % Azymuthal particle transport
+% heat flux
+Q_x = zeros(Nx,Ny,Nz,Ns3D);
+Q_y = zeros(Nx,Ny,Nz,Ns3D);
+
+% Epot averages
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_avgy = zeros(Nz,Ns3D); % Time evol. of the norm of phi
@@ -121,7 +127,13 @@ for it = 1:numel(Ts3D) % Loop over 2D aX_XYays
phi_avgx_avgy(iz,it) = mean(mean(squeeze(phi(:,:,iz,it))));
if(W_DENS)
- Gamma_x(:,:,iz,it) = dens_i(:,:,iz,it).*dyphi(:,:,iz,it);
+ Gamma_x(:,:,iz,it) = -dens_i(:,:,iz,it).*dyphi(:,:,iz,it);
+ Gamma_y(:,:,iz,it) = dens_i(:,:,iz,it).*dxphi(:,:,iz,it);
+ end
+
+ if(W_TEMP)
+ Q_x(:,:,iz,it) = -temp_e(:,:,iz,it).*dyphi(:,:,iz,it);
+ Q_y(:,:,iz,it) = temp_i(:,:,iz,it).*dxphi(:,:,iz,it);
end
shear_maxx_maxy(iz,it) = max( max(squeeze(-(dx2phi(:,:,iz,it)))));
diff --git a/matlab/setup.m b/matlab/setup.m
index f753c54dbd04ed5318f2e14ff32e372352f0297c..168c965ac96ea7962f3f1ea14ab61e2bc9a27e0a 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -55,9 +55,10 @@ if (abs(CO) == 2) %Sugama operator
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_NFLR_8.h5'''];
- INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_10_J_5_N_20_kpm_8.0_NFLR_0.h5'''];
-% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_k2.h5'''];
+% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5'''];
+ INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5'''];
+% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5'''];
+% INITIAL.mat_file = ['''../../../iCa/gk_coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5'''];
elseif (CO == -1) % DGDK
disp('Warning, DGDK not debugged')
end
@@ -71,7 +72,7 @@ switch abs(CO)
case 4; CONAME = 'FC';
otherwise; CONAME ='UK';
end
-if (CO < 0); CONAME = [CONAME,'DK'];
+if (CO <= 0); CONAME = [CONAME,'DK'];
else; CONAME = [CONAME,'GK'];
end
if (CLOS == 0); CLOSNAME = 'Trunc.';
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 3061a6179193704c5fc4883b26d8606571110667..0bdc0cff0e881659b351255674e601e521fcd955 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,29 +1,28 @@
addpath(genpath('../matlab')) % ... add
addpath(genpath('../matlab/plots')) % ... add
-%% Directory of the simulation
-if 1% Local results
-outfile ='';
-outfile ='';
-outfile ='';
outfile ='';
+%% Directory of the simulation
+if 0% Local results
outfile ='';
outfile ='';
-outfile ='HD_study/200x100_L_200_P_2_J_1_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_0e+00';
-% outfile ='simulation_B/cw_DGGK'; % to analyse
-% outfile ='simulation_B/cw_SGGK_mu_1e-1';
-% outfile ='simulation_B/cw_DGGK';
+% outfile ='artificial_ZF_freeze/sim_A';
+% outfile ='simulation_B/cw_FCGK_kp_3.0';
+% outfile ='nonlin_FCGK/150x75_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00';
+% outfile ='nonlin_PAGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_PAGK_mu_0e+00';
+% outfile ='nonlin_FCGK/100x50_L_200_P_4_J_2_eta_0.6_nu_1e-01_FCGK_mu_0e+00';
+outfile ='simulation_A';
+% outfile ='simulation_B/cw_SGGK_like_species';
+% outfile ='simulation_A/CO_damping_SGGK';
% outfile ='simulation_A/cw_DGGK_eta_0.5';
-% outfile ='simulation_A/LBDK_damping_150x75_L_100_P_4_J_2_eta_0.6_nu_1e-01_SGGK_mu_0e+00';
+% outfile ='simulation_B/cw_DGGK';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',BASIC.MISCDIR]);
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
system(CMD);
else% Marconi results
-outfile ='';
-outfile ='';
-outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt';
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
end
@@ -49,7 +48,7 @@ end
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1.2e4; TAVG_1 = 1.3e4; % Averaging times duration
+TAVG_0 = 1000; TAVG_1 = 2000; % Averaging times duration
plot_radial_transport_and_shear
end
@@ -63,37 +62,39 @@ end
if 0
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
% tstart = 0.8*Ts3D(end); tend = Ts3D(end); % Time window
-tstart = 14000; tend = 15000;
+tstart = 2500; tend = 2500;
% tstart = 10000; tend = 12000;
% Chose the field to plot
% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
-FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
+FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:76,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
-LOGSCALE = 0; TRENDS = 0; NORMALIZED = 0;
+LOGSCALE = 1; TRENDS = 1; NORMALIZED = 0;
plot_kperp_spectrum
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-t0 =0; iz = 1; ix = 1; iy = 1;
-skip_ =1; DELAY = 2e-3*skip_;
+t0 =000; iz = 1; ix = 1; iy = 1;
+skip_ =4; DELAY = 2e-3*skip_;
[~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
[~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
-INTERP = 0; T = Ts3D; FRAMES = FRAMES_3D;
+T = Ts3D; FRAMES = FRAMES_3D;
+INTERP = 0; NORMALIZED = 0; CONST_CMAP = 0;% Gif options
% Field to plot
% FIELD = dens_e; NAME = 'ne'; FIELDNAME = 'n_e';
% FIELD = dens_i; NAME = 'ni'; FIELDNAME = 'n_i';
% FIELD = dens_e-Z_n_e; NAME = 'ne_NZ';FIELDNAME = 'n_e^{NZ}';
-% FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
+FIELD = dens_i-Z_n_i; NAME = 'ni_NZ';FIELDNAME = 'n_i^{NZ}';
% FIELD = temp_e; NAME = 'Te'; FIELDNAME = 'n_i';
% FIELD = temp_i; NAME = 'Ti'; FIELDNAME = 'n_i';
% FIELD = temp_e-Z_T_e; NAME = 'Te_NZ';FIELDNAME = 'T_e^{NZ}';
% FIELD = temp_i-Z_T_i; NAME = 'Ti_NZ';FIELDNAME = 'T_i^{NZ}';
% FIELD = ne00; NAME = 'ne00'; FIELDNAME = 'n_e^{00}';
% FIELD = ni00; NAME = 'ni00'; FIELDNAME = 'n_i^{00}';
-FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = phi; NAME = 'phi'; FIELDNAME = '\phi';
+% FIELD = Z_phi; NAME = 'Zphi'; FIELDNAME = '\phi_Z';
% FIELD = Gamma_x; NAME = 'Gamma_x'; FIELDNAME = '\Gamma_x';
% Sliced
@@ -128,13 +129,14 @@ if 0
% FIELD = dens_i-Z_n_i; FNAME = 'ni_NZ'; FIELDLTX = 'n_i^{NZ}';
% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = 'T_i';
% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = 'T_e';
-FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
+% FIELD = phi; FNAME = 'phi'; FIELDLTX = '\phi';
% FIELD = Z_phi-phi; FNAME = 'phi_NZ'; FIELDLTX = '\phi^{NZ}';
+FIELD = Z_phi; FNAME = 'phi_Z'; FIELDLTX = '\phi^{Z}';
% FIELD = Gamma_x; FNAME = 'Gamma_x'; FIELDLTX = '\Gamma_x';
% FIELD = dens_e-Z_n_e-(Z_phi-phi); FNAME = 'Non_adiab_part'; FIELDLTX = 'n_e^{NZ}-\phi^{NZ}';
% Chose when to plot it
-tf = 128:1:133;
+tf = [9800 10000];
% Sliced
ix = 1; iy = 1; iz = 1;
@@ -157,7 +159,7 @@ plt_2 = @(x) x;%./max(max(x));
subplot(1,numel(tf),i_)
DATA = plt_2(squeeze(plt(FIELD,it)));
pclr = pcolor((X),(Y),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
- colormap(bluewhitered); caxis([-30,30]);
+ colormap(bluewhitered); %caxis([-30,30]);
xlabel(latexize(XNAME)); ylabel(latexize(YNAME));set(gca,'ytick',[]);
title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
end
@@ -169,14 +171,15 @@ if 0
%% Photomaton : k space
% Chose the field to plot
-FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
+% FIELD = Ni00; FNAME = 'Ni00'; FIELDLTX = 'N_i^{00}';
% FIELD = Ne00; FNAME = 'Ne00'; FIELDLTX = 'N_e^{00}'
-% FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
+FIELD = PHI; FNAME = 'PHI'; FIELDLTX = '\tilde\phi';
% FIELD_ = fft2(Gamma_x); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'Gamma_x'; FIELDLTX = '\tilde\Gamma_x';
% FIELD_ = fft2(dens_e); FIELD = FIELD_(1:Nx/2+1,:,:,:); FNAME = 'FFT_Dens_e'; FIELDLTX = '\tilde n_e';
% Chose when to plot it
-tf = 14000:50:14200;
+tf = 1000:500:3000;
+% tf = 8000;
% Sliced
ix = 1; iy = 1; iz = 1;
@@ -193,7 +196,7 @@ plt_2 = @(x) (fftshift(x,2));
subplot(1,numel(tf),i_)
DATA = plt_2(squeeze(plt(FIELD,it)));
pclr = pcolor(fftshift(X,2),fftshift(Y,2),DATA); set(pclr, 'edgecolor','none');pbaspect([0.5 1 1])
-% colormap(bluewhitered); caxis([-1,1]);
+ caxis([0 1]*5e3);
xlabel(latexize(XNAME)); ylabel(latexize(YNAME));
if(i_ > 1); set(gca,'ytick',[]); end;
title(sprintf('$t c_s/R=%.0f$',Ts3D(it)));
@@ -204,45 +207,104 @@ end
if 0
%%
-TAVG_0 = 10000; TAVG_1 = 15000; % Averaging times duration
+TAVG_0 = 1000; TAVG_1 = 5000; % Averaging times duration
ZF_fourier_analysis
end
-if 1
+
+if 0
%%
plot_param_evol
end
if 0
-%%
+%% Radial shear profile
+tf = 3000+[900:20:1100];
+ymax = 0;
figure
-plot(Ts3D,shear_maxx_avgy);
-
+for i_ = 1:numel(tf)
+[~,it] = min(abs(Ts3D-tf(i_)));
+data = squeeze((mean(dxphi(:,:,1,it),2)));
+plot(x,data,'Displayname',sprintf('$t c_s/R=%.0f$',Ts3D(it))); hold on;
+ymax = max([ymax abs(min(data)) abs(max(data))]);
+end
+xlim([min(x), max(x)]); ylim(1.2*[-1 1]*ymax);
+xlabel('$x/\rho_s$'); ylabel('$v_{E\times B,x}$'); grid on
end
-if 0
+if 1
%% zonal vs nonzonal energies for phi(t)
+trange = 200:Ns3D;
Ephi_Z = zeros(1,Ns3D);
-Ephi_NZ = zeros(1,Ns3D);
+Ephi_NZ_kgt0 = zeros(1,Ns3D);
+Ephi_NZ_kgt1 = zeros(1,Ns3D);
+Ephi_NZ_kgt2 = zeros(1,Ns3D);
+high_k_phi = zeros(1,Ns3D);
for it = 1:numel(Ts3D)
- Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
- Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
+% Ephi_NZ(it) = sum(sum(((KY~=0).*abs(PHI(:,:,1,it)).^2)));
+% Ephi_Z(it) = sum(sum(((KY==0).*abs(PHI(:,:,1,it)).^2)));
+ [amp,ikzf] = max(abs((kx~=0).*PHI(:,1,1,it)));
+% Ephi_NZ(it) = sum(sum(((KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+ Ephi_NZ_kgt0(it) = sum(sum(((sqrt(KX.^2+KY.^2)>0.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+ Ephi_NZ_kgt1(it) = sum(sum(((sqrt(KX.^2+KY.^2)>1.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+ Ephi_NZ_kgt2(it) = sum(sum(((sqrt(KX.^2+KY.^2)>2.0).*(KX~=0).*(KY~=0).*(KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)));
+% Ephi_Z(it) = kx(ikzf)^2*abs(PHI(ikzf,1,1,it)).^2;
+ Ephi_Z(it) = sum(sum(((KX~=0).*(KY==0).*(KX.^2).*abs(PHI(:,:,1,it)).^2)));
+% Ephi_NZ(it) = sum(sum(((KX.^2+KY.^2).*abs(PHI(:,:,1,it)).^2)))-Ephi_Z(it);
+ high_k_phi(it) = abs(PHI(18,18,1,it)).^2; % kperp sqrt(2)
+% high_k_phi(it) = abs(PHI(40,40,1,it)).^2;% kperp 3.5
end
pltx = @(x) x-x(1);
-plty = @(x) x./max(x);
+plty = @(x) x./max(squeeze(x));
fig = figure; FIGNAME = ['ZF_turb_energy_vs_time_',PARAMS];
set(gcf, 'Position', [100, 100, 1400, 500])
-subplot(131)
- semilogy(pltx(Ts3D),plty(Ephi_Z),'DisplayName',['$|\phi_k|^2$ ',CONAME]); hold on;
- title('Zonal Energy'); legend('show')
+subplot(121)
+% yyaxis left
+ semilogy(pltx(Ts3D(trange)),plty(Ephi_Z(trange)),'DisplayName',['Zonal, ',CONAME]); hold on;
+% yyaxis right
+ semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt0(trange)),'DisplayName',['NZ, $k_p>0$, ',CONAME]);
+ semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt1(trange)),'DisplayName',['NZ, $k_p>1$, ',CONAME]);
+ semilogy(pltx(Ts3D(trange)),plty(Ephi_NZ_kgt2(trange)),'DisplayName',['NZ, $k_p>2$, ',CONAME]);
+% semilogy(pltx(Ts0D),plty(PGAMMA_RI),'DisplayName',['$\Gamma_x$, ',CONAME]);
+ title('Energy'); legend('show')
xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
-subplot(132)
- semilogy(pltx(Ts3D),plty(Ephi_NZ));
- title('Non Zonal Energy'); legend(CONAME)
- xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
-
-subplot(133)
- semilogy(pltx(Ts0D),plty(abs(PGAMMA_RI)*SCALE));
- title('Radial particle flux'); legend(CONAME)
- xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
+subplot(122)
+ plot(plty(Ephi_Z(trange)),plty(Ephi_NZ_kgt0(trange)));
+ title('Phase space'); legend(CONAME)
+ xlabel('$E_Z$'); ylabel('$E_{NZ}$'); grid on;% xlim([0 500]);
+%
+% subplot(133)
+% % semilogy(pltx(Ts0D),plty(abs(PGAMMA_RI)*SCALE));
+% for ik = [10 20 30]
+% semilogy(pltx(Ts3D(trange)),plty(abs(squeeze(PHI(ik,ik,1,trange))).^2),'DisplayName',[CONAME,', kp=',num2str(sqrt(kx(ik)^2+ky(ik)^2))]); hold on
+% end
+% title('High kperp damping'); legend('show');
+% xlabel('$t c_s/R$'); grid on;% xlim([0 500]);
+end
+
+if 0
+%% Conservation laws
+Nxmax = Nx/2;
+Nymax = Ny/2;
+mflux_x_i = squeeze(sum((Gamma_x( 1,1:Nxmax,:)+Gamma_x( 1,2:Nxmax+1,:))/2,2)./sum(Gamma_x( 1,1:Nxmax,:)));
+mflux_x_o = squeeze(sum((Gamma_x( Nxmax,1:Nxmax,:)+Gamma_x( Nxmax,2:Nxmax+1,:))/2,2)./sum(Gamma_x( Nxmax,1:Nxmax,:)));
+mflux_y_i = squeeze(sum((Gamma_y(1:Nxmax, 1,:)+Gamma_y(2:Nxmax+1, 1,:))/2,1)./sum(Gamma_y(1:Nxmax, 1,:)));
+mflux_y_o = squeeze(sum((Gamma_y(1:Nxmax, Nymax,:)+Gamma_y(2:Nxmax+1, Nymax,:))/2,1)./sum(Gamma_y(1:Nxmax, Nymax,:)));
+
+mass_cons = mflux_x_i - mflux_x_o + mflux_y_i - mflux_y_o;
+%%
+figure
+plt = @(x) squeeze(mean(mean(x(:,:,1,:),1),2));
+subplot(211)
+ plot(Ts3D,plt(dens_e+dens_i),'DisplayName','$\delta n_e + \delta n_i$'); hold on;
+ plot(Ts3D,plt(ne00+ni00),'DisplayName','$\delta n_e^{00} + \delta n_i^{00}$'); hold on;
+ plot(Ts3D,plt(temp_e+temp_i),'DisplayName','$\delta T_e + \delta T_i$'); hold on;
+ legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)])
+subplot(212);
+ plot(Ts3D,mass_cons*(2*pi/Nx/Ny)^2,'DisplayName','in-out'); hold on
+% plot(Ts3D,squeeze(mflux_x_i),'DisplayName','Flux i x');
+% plot(Ts3D,squeeze(mflux_x_o),'DisplayName','Flux o x');
+% plot(Ts3D,squeeze(mflux_y_i),'DisplayName','Flux i y');
+% plot(Ts3D,squeeze(mflux_y_o),'DisplayName','Flux o y');
+ legend('show'); grid on; xlim([Ts3D(1) Ts3D(end)]); %ylim([-0.1, 2]*mean(mflux_x_i))
end
\ No newline at end of file
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index fc4fa8b74a005be18aa32f19d6465cb8d6ab8e1f..944063f38c24f180a09d8b755907689bda0dff84 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,4 +1,5 @@
%% Paste the list of continue_run calls
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
%% Functions to modify preexisting fort.90 input file and launch on marconi
@@ -44,19 +45,21 @@ function [] = continue_run(outfilename)
line = A{39};
line = line(end-2:end);
if(line(1) == '='); line = line(end); end;
- J2L = str2num(line) + 1;
+ J2L = str2num(line)+1;
end
% Change job 2 load in fort.90
A{39} = [' job2load = ',num2str(J2L)];
disp(A{39})
% Change time step
- A{3} = [' dt = 0.01'];
+ line_= A{3};
+ dt_old = str2num(line_(13:end));
+ A{3} = [' dt = ',num2str(dt_old)];
% Increase endtime
A{4} = [' tmax = 20000'];
% Change collision operator
line_= A{43};
CO_old = str2num(line_(13:end));
- A{43} = [' CO = ',num2str(1)];
+ A{43} = [' CO = ',num2str(2)];
% Put non linear term back
A{45} = [' NL_CLOS = -1'];
% change HD
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 58b08996205a97d8b824f2147f260fe00a917692..547bb0c58d65809259c17223bae679ffbe415790 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,4 +1,4 @@
-for CO = [2]
+for CO = [1]
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
@@ -12,7 +12,7 @@ TAU = 1.0; % e/i temperature ratio
ETAB = 1.0;
ETAN = 1/0.6; % Density gradient
ETAT = 0.0; % Temperature gradient
-NU_HYP = 3.0; % Hyperdiffusivity coefficient
+NU_HYP = 0.0; % Hyperdiffusivity coefficient
LAMBDAD = 0.0;
NOISE0 = 1.0e-5; % Init noise amplitude
BCKGD0 = 0.0; % Init background
@@ -24,11 +24,11 @@ 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 = 400; % Maximal time unit
-DT = 8e-3; % Time step
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/50; % Sampling per time unit for 5D arrays
+SPS3D = 2; % Sampling per time unit for 2D arrays
+SPS5D = 2; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 00;
@@ -36,13 +36,14 @@ JOB2LOAD= 00;
% SIMID = 'v3.6_kobayashi_lin'; % Name of the simulation
% SIMID = 'v3.2_CO_damping'; % Name of the simulation
% SIMID = 'CO_Patchwork_damping'; % Name of the simulation
-SIMID = 'v3.4_entropy_mode_linear'; % Name of the simulation
+SIMID = 'test_GF_closure'; % Name of the simulation
+% SIMID = 'v3.2_entropy_mode_linear'; % 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 = 1;
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)
+CLOS = 1; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
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
@@ -75,7 +76,7 @@ if 1
%% Parameter scan over PJ
% PA = [2 4];
% JA = [1 2];
-PA = [4];
+PA = [5];
JA = [2];
DTA= DT*ones(size(JA));%./sqrt(JA);
% DTA= DT;
@@ -86,8 +87,8 @@ param_name = 'PJ';
gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
gamma_Nipj = zeros(Nparam,floor(N/2)+1);
gamma_phi = zeros(Nparam,floor(N/2)+1);
-Ni00_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
- PHI_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
+% Ni00_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
+% PHI_ST = zeros(Nparam,floor(N/2)+1,TMAX/SPS3D);
for i = 1:Nparam
% Change scan parameter
PMAXE = PA(i); PMAXI = PA(i);
@@ -96,7 +97,8 @@ for i = 1:Nparam
setup
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3.6 1 6; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz; cd ../../../wk'])
end
% Load and process results
%%
@@ -135,7 +137,7 @@ plt = @(x) x;
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(1:numel(kx))),plt(gamma_Ni00(i,1:end)),...
+ plot(plt(SCALE*kx(1:numel(kx))),plt(gamma_Ni00(i,1:end)),...
'Color',clr,...
'LineStyle',linestyle{1},'Marker','^',...
...% 'DisplayName',['$\eta=',num2str(ETAB/ETAN),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
@@ -158,4 +160,11 @@ if 0
% pclr = pcolor(XT,YT,squeeze(abs(Ni00_ST(1,:,:)))); set(pclr, 'edgecolor','none'); colorbar;
semilogy(Ts3D(1:TMAX/SPS3D),squeeze(abs(PHI_ST(1,50:5:100,:))));
end
+if 0
+%% Trajectories of some modes
+figure;
+for i = 1:10:N/2+1
+ semilogy(Ts3D,squeeze(abs(Ne00(i,1,1,:))),'DisplayName',['k=',num2str(kx(i))]); hold on;
+end
+end
end
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index 87fe34405d3dddf761422a0cc24658d362ea2b49..0738cad9421e410b82e9d123eb6b13a1d2869ec2 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,6 +1,4 @@
addpath(genpath('../matlab')) % ... add
%% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_DGGK_mu_0e+00/out.txt')
load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_DGGK_mu_0e+00/out.txt')
load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_B/300x150_L_120_P_8_J_4_eta_0.6_nu_5e-01_SGGK_mu_0e+00/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 542b50918348241e0f9fcac834699ec89cabd6ef..7de7329dfa4c3dfe2fe6bfb26c699c638f2832fc 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -6,10 +6,10 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
ETAN = 1.0/0.6; % Density gradient drive (R/Ln)
-NU_HYP = 1.0;
+NU_HYP = 0.0;
%% GRID AND GEOMETRY PARAMETERS
N = 150; % Frequency gridpoints (Nkx = N/2)
-L = 100; % Size of the squared frequency domain
+L = 200; % Size of the squared frequency domain
Nz = 1; % number of perpendicular planes (parallel grid)
q0 = 1.0; % safety factor
shear = 0.0; % magnetic shear
@@ -17,25 +17,25 @@ eps = 0.0; % inverse aspect ratio
P = 4;
J = 2;
%% TIME PARAMETERS
-TMAX = 1000; % Maximal time unit
-DT = 1e-3; % Time step
+TMAX = 100; % Maximal time unit
+DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS3D = 1/2; % Sampling per time unit for 3D arrays
+SPS3D = 1; % Sampling per time unit for 3D arrays
SPS5D = 1/20; % 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;
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; 4 : Coulomb; +/- for GK/DK)
+CO = 4;
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 = 'HD_study'; % Name of the simulation
+NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
+SIMID = 'nonlin_FCGK'; % Name of the simulation
% SIMID = 'test_3D'; % Name of the simulation
% SIMID = ['v3.0_P_',num2str(P),'_J_',num2str(J)]; % Name of the simulation
-NON_LIN = -1; % activate non-linearity (is cancelled if KXEQ0 = 1)
+NON_LIN = 1; % activate non-linearity (is cancelled if KXEQ0 = 1)
% INIT options
INIT_ZF = 0; ZF_AMP = 0.0;
INIT_BLOB = 0; WIPE_TURB = 0; WIPE_ZF = 0;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 1620081fad564ff2af6cd3a96dee5554c77e738f..577e6b5aebf065d757998c836d5bcaf1a4afd909 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -2,7 +2,7 @@ clear all;
addpath(genpath('../matlab')) % ... add
SUBMIT = 1; % To submit the job automatically
% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_3.2';
+ EXECNAME = 'helaz_3.8';
for ETAN = [1/0.6]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index 60c8e489fa1f1050f0bd5949b6bedbe7836d045c..bade9a4d4e2059b2db6671e0941f769b323eda75 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -93,8 +93,12 @@ subplot(224)
%% Single eigenvalue analysis
% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_bjf/scanfiles_00000/self.0.h5';
-mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_6_3_NFLR_8_kp_4_new_Tljpmf_fort_bjf/scanfiles_00000/self.0.h5';
-matidx = 74;
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_4_2_NFLR_2_kp_4/scanfiles_00000/self.0.h5';
+mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0/scanfiles_00049/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_16_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
+% mat_file_name = '/home/ahoffman/cosolver/trunc_gk.coulomb_NFLR_4_P_4_J_2_N_50_kpm_4.0/scanfiles_00049/self.0.h5';
+matidx = 01;
matidx = sprintf('%5.5i',matidx);disp(matidx);
@@ -105,19 +109,21 @@ gmax = max(real(eig(MAT)));
wmax = max(imag(eig(MAT)));
figure
+subplot(121)
imagesc((MAT)); colormap(bluewhitered)
title(['$\gamma_{max}=',num2str(gmax),'$'])
+subplot(122)
+plot(real(eig(MAT)),imag(eig(MAT)),'x')
%% Eigenvalue spectrum 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_10_J_5_N_20_kpm_8.0_NFLR_0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_2_J_1_N_20_kpm_8.0_NFLR_5.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0_NFLR_8.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_75_kpm_6.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
+ '/home/ahoffman/cosolver/gk.coulomb_NFLR_16_P_4_J_2_N_50_kpm_4.0.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 10 5 NFLR 0', 'FC 2 1 NFLR 5', 'FC 6 3 NFLR 0', 'FC 6 3 NFLR 8'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 16'};
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
@@ -130,15 +136,15 @@ for j_ = 1:numel(mfns)
gmax(idx_+1) = max(real(eig(MAT)));
wmax(idx_+1) = max(imag(eig(MAT)));
end
-% subplot(121)
+ subplot(121)
plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
-% subplot(122)
-% plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
+ subplot(122)
+ plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
end
-% subplot(121)
+ 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)')
+ subplot(122)
+legend('show'); grid on;
+xlabel('$k_\perp$'); ylabel('$\omega_{max}$ from Eig(iCa)')
end