From f225dc1f7a28d2e8318afc7289dd30e99a6cdc89 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 15 Mar 2022 08:25:58 +0100
Subject: [PATCH] scripts updates
---
matlab/COSOlver_matrices_analysis.m | 72 -------------------
matlab/assemble_cosolver_matrices.m | 71 ++++++++++++++++++
matlab/extract_fig_data.m | 18 ++---
matlab/load/load_params.m | 4 +-
matlab/load_write_mat.m | 15 ++++
matlab/plot/plot_cosol_mat.m | 21 +++---
matlab/plot/plot_param_evol.m | 8 +++
.../plot_radial_transport_and_spacetime.m | 32 +++++----
matlab/plot/statistical_transport_averaging.m | 8 +--
matlab/setup.m | 5 +-
wk/analysis_3D.m | 34 ++++-----
wk/analysis_header.m | 13 +++-
wk/benchmark_HeLaZ_gene_transport_zpinch.m | 34 +++++----
wk/linear_1D_entropy_mode.m | 18 ++---
wk/mode_growth_meter.m | 2 +-
wk/quick_gene_load.m | 5 +-
16 files changed, 199 insertions(+), 161 deletions(-)
delete mode 100644 matlab/COSOlver_matrices_analysis.m
create mode 100644 matlab/assemble_cosolver_matrices.m
create mode 100644 matlab/load_write_mat.m
diff --git a/matlab/COSOlver_matrices_analysis.m b/matlab/COSOlver_matrices_analysis.m
deleted file mode 100644
index 083eb5ea..00000000
--- a/matlab/COSOlver_matrices_analysis.m
+++ /dev/null
@@ -1,72 +0,0 @@
-addpath(genpath('../matlab')) % ... add
-%% Grid configuration
-N = 10; % Frequency gridpoints (Nkx = N/2)
-L = 120; % Size of the squared frequency domain
-dk = 2*pi/L;
-kmax = N/2*dk;
-kx = dk*(0:N/2);
-ky = dk*(0:N/2);
-[ky, kx]= meshgrid(ky,kx);
-KPERP = sqrt(kx.^2 + ky.^2);
-kperp = reshape(KPERP,[1,numel(kx)^2]);
-kperp = uniquetol(kperp,1e-14);
-Nperp = numel(kperp);
-COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
-COSOLVER.pmaxe = 10;
-COSOLVER.jmaxe = 5;
-COSOLVER.pmaxi = 10;
-COSOLVER.jmaxi = 5;
-
-COSOLVER.neFLR = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation
-COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2));
-
-COSOLVER.neFLRs = 0; % ... only for GK abel
-COSOLVER.npeFLR = 0; % ... only for GK abel
-COSOLVER.niFLRs = 0; % ... only for GK abel
-COSOLVER.npiFLR = 0; % ... only for GK abel
-
-COSOLVER.gk = 1;
-COSOLVER.co = 3;
-if 0
- %% plot the kperp distribution
- figure
- plot(kperp)
-end
-%% Load the matrices
-C_self_i = zeros((COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),(COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),Nperp);
-
-for n_ = 1:Nperp
- COSOLVER.kperp = kperp(n_);
- k_string = sprintf('%0.4f',kperp(n_));
- mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
- '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
- '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
- '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
- '_JE_12',...
- '_NFLR_',num2str(COSOLVER.neFLR),...
- '_kperp_',k_string,'.h5'];
-
- tmp = h5read(mat_file_name,'/Caapj/Ceepj');
- C_self_i(:,:,n_) = tmp;
-end
-
-%% Post processing
-dC_self_i = diff(C_self_i,1,3);
-gvar_dC = squeeze(sum(sum(abs(dC_self_i),2),1));
-%% Plots
-%% all coeffs evolution
-figure
-for ip_ = 1:COSOLVER.pmaxi+1
- for ij_ = 1:COSOLVER.jmaxi+1
- plot(kperp,squeeze(C_self_i(ip_,ij_,:)),'o'); hold on;
- end
-end
-%% global matrix variation
-figure;
-kperpmid = 0.5*(kperp(2:end)+kperp(1:end-1));
-plot(kperpmid,gvar_dC./diff(kperp)');
-
-
-
-
-
diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
new file mode 100644
index 00000000..2dad7e0b
--- /dev/null
+++ b/matlab/assemble_cosolver_matrices.m
@@ -0,0 +1,71 @@
+codir = '/home/ahoffman/cosolver/';
+matname= 'gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0';
+matdir = [codir,matname,'/'];
+kperp = load([matdir,'kperp.log']);
+Nkp = numel(kperp);
+outdir = '/home/ahoffman/HeLaZ/iCa/';
+outfile= [outdir,matname,'.h5'];
+
+if(exist(outfile)>0)
+ system(['rm ',outdir,matname,'.h5']);
+end
+
+Nmax = numel(kperp);
+for n_=1:Nmax
+ n_string = sprintf('%5.5d',n_-1); disp(n_string);
+ scandir = ['scanfiles_',n_string,'/'];
+ if(exist([matdir,scandir, 'ei.h5'])>0)
+ % Load matrices and write them in one h5 file
+
+ olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
+ newdname = ['/',n_string,olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ceipj/CeipjT'; infile = [matdir,scandir, 'ei.h5'];
+ newdname = ['/',n_string,olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ciepj/CiepjT'; infile = [matdir,scandir, 'ie.h5'];
+ newdname = ['/',n_string,olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Ciepj/CiepjF'; infile = [matdir,scandir, 'ie.h5'];
+ newdname = ['/',n_string,olddname];
+ load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Caapj/Ceepj'; infile = [matdir,scandir,'self.h5'];
+ newdname = ['/',n_string,olddname];
+ mat_ee = load_write_mat(infile,olddname,outfile,newdname);
+
+ olddname = '/Caapj/Ciipj'; infile = [matdir,scandir,'self.h5'];
+ newdname = ['/',n_string,olddname];
+ mat_ii = load_write_mat(infile,olddname,outfile,newdname);
+
+ % to verify symmetry
+ verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii))));
+ if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
+ % to verify negative EV
+ verif = max(real(eig(mat_ee))) + max(real(eig(mat_ii)));
+ if(verif>0) disp(['Warning, positive EV at ', n_string]); end;
+ else
+ break
+ end
+end
+olddname = '/Ceipj/CeipjF'; infile = [matdir,scandir, 'ei.h5'];
+Pmaxe = h5readatt(infile,olddname,'Pmaxe');
+Jmaxe = h5readatt(infile,olddname,'Jmaxe');
+Pmaxi = h5readatt(infile,olddname,'Pmaxi');
+Jmaxi = h5readatt(infile,olddname,'Jmaxi');
+dims_e = [0 0];
+dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
+dims_i = [Pmaxi; Jmaxi];
+h5create(outfile,'/dims_e',numel(dims_e));
+h5write (outfile,'/dims_e',dims_e);
+h5create(outfile,'/dims_i',numel(dims_i));
+h5write (outfile,'/dims_i',dims_i);
+%
+h5create(outfile,'/coordkperp',numel(kperp(1:Nmax)));
+h5write (outfile,'/coordkperp',kperp(1:Nmax)');
+fid = H5F.open(outfile);
+
+H5F.close(fid);
\ No newline at end of file
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index e82fcf97..27f177e4 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -7,20 +7,20 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData];
end
-n0 = 100;
+n0 = 1;
figure;
mvm = @(x) movmean(x,1);
shift = X_(n0);
% shift = 0;
% plot(X_(n0:end),Y_(n0:end));
-plot(mvm(X_(n0:end)),mvm(Y_(n0:end))./max(Y_(200:end))); hold on;
+plot(mvm(X_(n0:end)-shift),mvm(Y_(n0:end))); hold on;
-
-n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
-avg_ = mean(Y_(n1:n2));
-std_ = std(Y_(n1:n2));
-title(['avg = ',num2str(avg_),' std = ', num2str(std_)])
-plot([X_(n1),X_(n2)],[1 1]*avg_,'--k')
+%
+% n1 = n0+1; n2 = min(n1 + 50000,numel(Y_));
+% avg_ = mean(Y_(n1:n2));
+% std_ = std(Y_(n1:n2));
+% title(['avg = ',num2str(avg_),' std = ', num2str(std_)])
+% plot([X_(n1),X_(n2)],[1 1]*avg_,'--k')
% ylim([0,avg_*3]);
sum(Y_(2:end)./X_(2:end).^(3/2))
@@ -35,4 +35,4 @@ if 0
end
- xlim([-3,3]); ylim([0,4.5]);
\ No newline at end of file
+% xlim([-3,3]); ylim([0,4.5]);
\ No newline at end of file
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 78020b7f..613a62bd 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -20,8 +20,8 @@ DATA.L = h5readatt(filename,'/data/input','Lx');
DATA.CLOS = h5readatt(filename,'/data/input','CLOS');
DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
-DATA.MUx = -99;%h5readatt(filename,'/data/input','mu_x');
-DATA.MUy = -99;%h5readatt(filename,'/data/input','mu_y');
+DATA.MUx = h5readatt(filename,'/data/input','mu_x');
+DATA.MUy = h5readatt(filename,'/data/input','mu_y');
% MU = str2num(filename(end-18:end-14)); %bad...
DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
diff --git a/matlab/load_write_mat.m b/matlab/load_write_mat.m
new file mode 100644
index 00000000..3315b172
--- /dev/null
+++ b/matlab/load_write_mat.m
@@ -0,0 +1,15 @@
+function mat_ = load_write_mat(infile,olddname,outfile,newdname)
+ mat_ = h5read (infile,olddname);
+
+ h5create (outfile,newdname,size(mat_));
+ h5write (outfile,newdname,mat_);
+
+ if(~strcmp(olddname(end-3:end-2),'ii'))
+ neFLR = h5readatt(infile,olddname,'neFLR');
+ h5writeatt(outfile,newdname,'neFLR',neFLR);
+ end
+ if(~strcmp(olddname(end-3:end-2),'ee'))
+ niFLR = h5readatt(infile,olddname,'niFLR');
+ h5writeatt(outfile,newdname,'niFLR',niFLR);
+ end
+end
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index faa3c152..6f576a1d 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -116,11 +116,12 @@ 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/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 16'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG A', 'Hacked SG B'};
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
@@ -150,14 +151,14 @@ end
%% Van Kampen plot
if 0
%%
-kperp= 3.9;
+kperp= 0.1;
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/cosolver/gk.coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/HeLaZ/wk/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 4 2 NFLR 12 extended'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'Hacked SG'};
grow = {};
puls = {};
for j_ = 1:numel(mfns)
@@ -177,4 +178,4 @@ for j_ = 1:numel(mfns)
end
legend('show'); grid on;
xlabel('$\omega$ from Eig(iCa)'); ylabel('$\gamma$ from Eig(iCa)')
-end
+end
\ No newline at end of file
diff --git a/matlab/plot/plot_param_evol.m b/matlab/plot/plot_param_evol.m
index b8707232..e1f13b2c 100644
--- a/matlab/plot/plot_param_evol.m
+++ b/matlab/plot/plot_param_evol.m
@@ -9,12 +9,20 @@ MUy_EVOL = data.MUy_EVOL;
K_N_EVOL = data.K_N_EVOL;
L_EVOL = data.L_EVOL;
+CO_LIST = {};
+for i_ = 1:numel(CO_EVOL)/6
+ CO_LIST{i_} = CO_EVOL(6*(i_-1)+1:6*i_);
+end
fig = figure;
hold on; set(gcf, 'Position', [100, 100, 400, 1500])
subplot(4,1,1);
title('Parameter evolution'); hold on;
yyaxis left
plot(TJOB_SE,NU_EVOL,'DisplayName','$\nu$');
+ for i_ = 1:numel(CO_EVOL)/12
+ k_ = 2*(i_-1)+1;
+ text(double(TJOB_SE(k_)),NU_EVOL(k_)*2.0,CO_LIST{k_});
+ end
yyaxis right
plot(TJOB_SE,DT_EVOL,'DisplayName','dt');
xlim([TJOB_SE(1) TJOB_SE(end)]);legend('show');grid on;
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index fa88f634..9a3234df 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -45,21 +45,6 @@ mvm = @(x) movmean(x,Nmvm);
grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
title([DATA.param_title,', $\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
- % plot
- subplot(312)
- it0 = 1; itend = Ns3D;
- trange = it0:itend;
- plt1 = @(x) x;%-x(1);
- plt2 = @(x) x./max((x(its3D:ite3D)));
-% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
- yyaxis left
- plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
- ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
- yyaxis right
- plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
- xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
- ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
- xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
@@ -108,4 +93,21 @@ mvm = @(x) movmean(x,Nmvm);
subplot(311)
plot(DATA.Ts3D,squeeze(mean(plt(f2plot),1)));
end
+%% Zonal vs NZonal energies
+ subplot(312)
+ it0 = 1; itend = Ns3D;
+ trange = it0:itend;
+ plt1 = @(x) x;%-x(1);
+ plt2 = @(x) x./max((x(its3D:ite3D)));
+ toplot = sum(squeeze(plt(f2plot)).^2,1); %ST from before
+% plty = @(x) x(500:end)./max(squeeze(x(500:end)));
+ yyaxis left
+% plot(plt1(DATA.Ts3D(trange)),plt2(E_Zmode_SK(trange)),'DisplayName','$k_{zf}^2|\phi_{kzf}|^2$');
+ plot(plt1(DATA.Ts3D(trange)),plt2(toplot(trange)),'DisplayName','Sum $A^2$');
+ ylim([-0.1, 1.5]); ylabel('$E_{Z}$')
+ yyaxis right
+ plot(plt1(DATA.Ts3D(trange)),plt2(E_NZmode_SK(trange)),'DisplayName','$(1+k^2)|\phi^k|^2$');
+ xlim([DATA.Ts3D(it0), DATA.Ts3D(itend)]);
+ ylim([-0.1, 1.5]); ylabel('$E_{NZ}$')
+ xlabel('$t c_s/R$'); grid on; set(gca,'xticklabel',[]);% xlim([0 500]);
end
\ No newline at end of file
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index 5a9d0653..67ed61da 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -5,9 +5,9 @@ Trange = options.T;
[~,it1] = min(abs(Trange(end)-data.Ts0D));
gamma = data.PGAMMA_RI(it0:it1)*scale;
dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
-if ~dt_const
- disp('DT not const on given interval');
-else
+% if ~dt_const
+% disp('DT not const on given interval');
+% else
Ntot = (it1-it0)+1;
@@ -29,6 +29,6 @@ else
% subplot(212)
% errorbar(N_seg,transp_seg_avg,transp_seg_std);
% xlabel('Averaging #points'); ylabel('transport value');
-end
+% end
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 5ec853a9..d1afb57a 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -51,10 +51,13 @@ COLL.mat_file = '''null''';
switch CO
case 'SG'
COLL.mat_file = '''../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5''';
+% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5''';
+% COLL.mat_file = '''../../../iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5''';
case 'LR'
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD'
- COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''';
+% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_75_kpm_6.0.h5''';
+ COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
end
% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 568a56d7..cb59d386 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -10,7 +10,7 @@ system(['mkdir -p ',MISCDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 10;
+JOBNUMMIN = 00; JOBNUMMAX = 20;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -21,7 +21,7 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 500; TAVG_1 = 10000; % Averaging times duration
+TAVG_0 = 1000; TAVG_1 = 11000; % Averaging times duration
% chose your field to plot in spacetime diag (uzf,szf,Gx)
fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'v_y',1);
save_figure(data,fig)
@@ -29,7 +29,7 @@ end
if 0
%% statistical transport averaging
-options.T = [5000 6000];
+options.T = [1500 2500];
fig = statistical_transport_averaging(data,options);
end
@@ -39,16 +39,16 @@ if 0
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
-% options.NAME = 'v_y';
+% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+% options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 1;
% options.TIME = data.Ts5D;
-options.TIME = 0:3:1000;
+options.TIME = 900:10:2000;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
end
@@ -56,19 +56,19 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 0;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
% options.NAME = '\phi';
-% options.NAME = 'v_y';
-options.NAME = 'n_e^{NZ}';
+options.NAME = 'v_y';
+% options.NAME = 'n_e^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [550:20:800];
+options.TIME = [000:200:1000];
data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -92,7 +92,7 @@ options.XPERP = linspace( 0,6,64);
% options.SPAR = vp';
% options.XPERP = mu';
options.Z = 1;
-options.T = 4000;
+options.T = 5500;
options.CTR = 1;
options.ONED = 0;
fig = plot_fa(data,options);
@@ -104,7 +104,7 @@ if 0
% options.TIME = 'avg';
options.TIME = 1000:4000;
options.P2J = 1;
-options.ST = 1;
+options.ST = 0;
options.NORMALIZED = 0;
fig = show_moments_spectrum(data,options);
save_figure(data,fig)
@@ -112,15 +112,15 @@ end
if 0
%% 1D spectrum
-options.TIME = 5000:10:5200;
+options.TIME = 5000:10:5050;
options.NORM = 1;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'n_i';
-options.NAME ='\Gamma_x';
+% options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 1;
options.OK = 0;
-options.COMPXY = 'max';
+options.COMPXY = 'sum';
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
@@ -147,7 +147,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 0.01:0.01:1.0;
-options.TIME = 5000:1:5050;
+options.TIME = 4900:1:5100;
options.NMA = 1;
options.NMODES = 20;
fig = mode_growth_meter(data,options);
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index e3f8e09f..4ac8664b 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -120,15 +120,18 @@ outfile ='';
% outfile = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5';
+% outfile = 'nu_0.1_transport_scan/LB_kn_2.0';
% outfile = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1';
% outfile = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5';
+% outfile = 'nu_0.1_transport_scan/DG_conv_kN_1.9';
% outfile = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/SG_10x5_conv_test';
% outfile = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5';
% outfile = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5';
+% outfile = 'nu_0.1_transport_scan/LD_kn_1.7_to_2.5';
% outfile = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0';
% outfile = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5';
@@ -139,10 +142,16 @@ outfile ='';
% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
% outfile = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1';
+% outfile = 'nu_0.1_transport_scan/large_box_kN_2.0_nu_0.1';
-% outfile = 'predator_prey/SG_Kn_1.7_nu_0.01';
+outfile = 'predator_prey_nu_scan/DG_Kn_1.7_nu_0.01';
-outfile = 'ZF_damping_DK/SG_4x2_nu_0.1';
+% outfile = 'ZF_damping_linear_nu_0_20x10_kn_1.6_GK/LR_4x2_nu_0.1';
+% outfile = 'ZF_damping_nu_0_20x10_kn_1.6_GK/HSG_4x2_nu_0.1';
+% outfile = 'ZF_damping_nu_0_5x3_kn_2.5_GK/LR_4x2_nu_0.1';
+% outfile = 'hacked_sugama/hacked_B_kn_1.6_200x32_L_120x60_nu_0.1';
+
+% outfile = 'shearless_cyclone/200x32x24_5x4_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_0';
% else% Marconi results
% outfile ='';
% outfile ='';
diff --git a/wk/benchmark_HeLaZ_gene_transport_zpinch.m b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
index 59d756eb..f0e5ddc6 100644
--- a/wk/benchmark_HeLaZ_gene_transport_zpinch.m
+++ b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
@@ -83,20 +83,14 @@ fig = figure; set(gcf,'Position',[250 250 600 300]);
% SGGK_200x32x11x06 = [1.5e-2 1.3e+0 0.0e+0];
% DGGK_200x32x11x06 = [2.4e-3 3.8e-1 0.0e+0];
-% upward scan
-KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
-LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.3e-1 1.3e+0 2.3e+0 3.6e+0 6.3e+0 9.4e+0 1.5e+1 2.3e+1 3.8e+1];
-LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.3e+0 3.5e+0 6.0e+0 4.8e+0 6.0e+0 2.9e+1];
-SGGK_200x32x05x03 = [0.0e-9 1.5e-2 1.0e-1 3.8e-1 8.4e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
-DGGK_200x32x05x03 = [2.0e-4 3.6e-3 1.3e-2 4.0e-2 1.2e-1 1.8e-1 2.0e-1 2.2e+0 6.6e-1 1.3e+1 3.2e+1];
-NOCO_200x32x05x03 = [0.0e+0 1.2e-2 2.6e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 3.9e+0 6.9e+0 1.2e+1 2.4e+1];
-% downward scan
-% KN = [KN 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
-% LDGK_200x32x05x03 = [LDGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-% LRGK_200x32x05x03 = [LRGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-% SGGK_200x32x05x03 = [SGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-% DGGK_200x32x05x03 = [DGGK_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-% NOCO_200x32x05x03 = [NOCO_200x32x05x03 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 3.9e+0 6.9e+0 1.2e+1 0.0e+0];
+KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5];
+LDGK_200x32x05x03 = [2.6e-2 2.6e-1 6.2e-1 1.2e+0 2.3e+0 3.6e+0 7.4e+0 1.1e+1 1.5e+1 2.3e+1 3.8e+1]; % validate with two box sizes (120 and 240)
+LRGK_200x32x05x03 = [0.0e+0 4.0e-1 7.2e-1 1.2e+0 1.7e+0 2.3e+0 4.0e+0 6.0e+0 1.1e+1 1.9e+1 2.9e+1];
+SGGK_200x32x05x03 = [0.0e-9 1.5e-2 9.5e-2 3.5e-1 8.0e-1 1.4e+0 2.5e+0 4.2e+0 8.0e+0 1.6e+1 3.5e+1];
+hacked_sugama = [0.0e+0 7.0e-1 0.0e+0 2.5e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+DGGK_200x32x05x03 = [2.0e-4 5.0e-3 9.2e-3 3.5e-2 1.3e-1 3.5e-1 6.0e-1 6.0e+0 1.1e+1 1.3e+1 3.2e+1]; % validate with two box sizes (120 and 240)
+LBGK_200x32x05x03 = [0.0e+0 0.0e-1 0.0e+0 0.0e+0 0.0e+0 2.8e-1 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0]; % validate with two box sizes (120 and 240)
+NOCO_200x32x05x03 = [0.0e+0 1.2e-2 4.0e-2 9.0e-2 2.3e-1 3.1e-1 7.0e-1 3.9e+0 6.9e+0 1.2e+1 2.4e+1];
% X = Ginf_GENE;
@@ -109,6 +103,7 @@ semilogy(X,plt(DGGK_200x32x05x03),'v','DisplayName','Dougherty','MarkerSize',msi
semilogy(X,plt(SGGK_200x32x05x03),'s','DisplayName', 'Sugama','MarkerSize',msize,'Color',clr_(1,:)); hold on;
semilogy(X,plt(LDGK_200x32x05x03),'d','DisplayName', 'Landau','MarkerSize',msize,'Color',clr_(5,:)); hold on;
semilogy(X,plt(LRGK_200x32x05x03),'^','DisplayName', 'Lorentz','MarkerSize',msize,'Color',clr_(3,:)); hold on;
+semilogy(X,plt(hacked_sugama),'s','DisplayName', 'hacked Sugama','MarkerSize',msize,'Color',clr_(6,:)); hold on;
semilogy(X,plt(NOCO_200x32x05x03),'*k','DisplayName', '$\nu = 0$','MarkerSize',msize); hold on;
X = [ 1.6 2.0 2.5];
@@ -122,11 +117,14 @@ xlim([ 1.55 2.55]);
% saveas(fig,'/home/ahoffman/Dropbox/Applications/Overleaf/Paper 1/figures/coll_transport_benchmark.eps')
%% Burst study Kn = 1.7
+Kn = 1.7;
NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 6.0e-2 7.0e-2 8.0e-2 9.0e-2 1.0e-1];
SGGK_200x32x05x03 = [1.0e-2 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 6.5e-2 7.5e-2 8.3e-2 8.8e-2 1.0e-1];
-SGGK_200x32x10x05 = [1.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-LRGK_200x32x05x03 = [0.0e+0 2.9e-2 0.0e+0 0.0e-2 0.0e-2 0.0e+0 0.0e+0];
+DGGK_200x32x05x03 = [1.0e-2 7.0e-3 7.3e-3 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+LDGK_200x32x05x03 = [1.0e-2 8.5e-2 1.4e-1 2.0e-1 0.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
% NOCO_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
-plot(NU,SGGK_200x32x05x03,'x'); hold on
-grid on; xlabel('$\nu$'); ylabel('$\Gamma^\infty_x$');
+plot(NU,SGGK_200x32x05x03/Kn,'s','Color',clr_(1,:)); hold on
+plot(NU,DGGK_200x32x05x03/Kn,'v','Color',clr_(2,:)); hold on
+plot(NU,LDGK_200x32x05x03/Kn,'d','Color',clr_(5,:)); hold on
+grid on; xlabel('$\nu$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
diff --git a/wk/linear_1D_entropy_mode.m b/wk/linear_1D_entropy_mode.m
index e6d1265e..1a74761d 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/wk/linear_1D_entropy_mode.m
@@ -8,12 +8,12 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 1.7; % Density gradient drive
+K_N = 1.6; % Density gradient drive
K_T = 0.25*K_N; % Temperature '''
K_E = 0.0; % Electrostat '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
%% GRID PARAMETERS
-NX = 64; % real space x-gridpoints
+NX = 40; % real space x-gridpoints
NY = 1; % '' y-gridpoints
LX = 120; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
@@ -23,8 +23,8 @@ SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 200; % Maximal time unit
+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
@@ -37,11 +37,11 @@ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
KIN_E = 1;
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'DG';
+CO = 'LD';
GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
-CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))
+CLOS = 0; % Closure model (0: =0 truncation, 1: gyrofluid closure (p+2j<=Pmax))s
NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
@@ -72,11 +72,11 @@ CURVB = 1.0;
if 1
% Parameter scan over PJ
-% PA = [4 6 8 10];
-% JA = [2 3 4 5];
+PA = [4];
+JA = [2];
Nparam = numel(PA);
% Parameter scan over KN
-PA = [4]; JA = [2];
+% PA = [4]; JA = [2];
% PMAXE = PA(1); PMAXI = PA(1);
% JMAXE = JA(1); JMAXI = JA(1);
% KNA = 1.5:0.05:2.5;
diff --git a/wk/mode_growth_meter.m b/wk/mode_growth_meter.m
index 5ff7721e..95af8e71 100644
--- a/wk/mode_growth_meter.m
+++ b/wk/mode_growth_meter.m
@@ -69,7 +69,7 @@ end
subplot(2,3,1+3*(i-1))
[YY,XX] = meshgrid(t,fftshift(k,1));
pclr = pcolor(XX,YY,abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
-% pclr = imagesc(t,fftshift(k,1),abs(plt(fftshift(DATA.PHI,MODES_SELECTOR),1:numel(k))));
+ set(gca,'YDir','normal')
% xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
title(MODESTR)
diff --git a/wk/quick_gene_load.m b/wk/quick_gene_load.m
index 7fd7a569..3eddddf2 100644
--- a/wk/quick_gene_load.m
+++ b/wk/quick_gene_load.m
@@ -24,6 +24,9 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_1.5_KT_0.325_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.4_KT_0.35_nuSG_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSG_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.6_KT_0.4_nuSGDK_0.0047_64x32.txt';
+fname ='GENE_LIN_Kn_1.6_KT_0.4_nuLDDK_0.0047_64x32.txt';
+% fname ='GENE_LIN_Kn_1.8_KT_0.4_nuLDDK_0.0047_64x32.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_1.8_KT_0.45_nuSG_0.0235_64x32.txt';
% fname ='GENE_LIN_Kn_1.9_KT_0.475_nuSG_0.047_64x32.txt';
@@ -36,7 +39,7 @@ path = '/home/ahoffman/gene/linear_zpinch_results/';
% fname ='GENE_LIN_Kn_2.0_KT_0.5_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.0_KT_0.5_nuSGDK_0.0235_32x16.txt';
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
-fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
+% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
data_ = load([path,fname]);
figure
--
GitLab