From 6571d244780a1ded603b99ecfb656a22e4b5a62f Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 2 Dec 2022 13:49:53 +0100
Subject: [PATCH] save the scripts
---
.../compute/compute_3D_zpinch_growth_rate.m | 10 +++-
matlab/compute/mode_growth_meter.m | 53 ++++++++++-------
matlab/dbg_zBC_map.m | 38 +++++++++----
matlab/evolve_tracers.m | 8 +--
matlab/extract_fig_data.m | 2 +-
matlab/load/load_gene_data.m | 2 +-
.../plot_radial_transport_and_spacetime.m | 3 +-
matlab/plot/spectrum_1D.m | 2 +-
matlab/process_field.m | 2 +-
matlab/setup.m | 2 +
matlab/write_fort90.m | 1 +
wk/analysis_gbms.m | 4 +-
wk/analysis_gene.m | 14 ++++-
wk/analysis_gyacomo.m | 57 +++++++++++--------
wk/header_2DZP_results.m | 9 ++-
wk/header_3D_results.m | 27 ++++++---
wk/lin_ITG.m | 32 ++++++-----
17 files changed, 168 insertions(+), 98 deletions(-)
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
index ce129abd..521e5a1f 100644
--- a/matlab/compute/compute_3D_zpinch_growth_rate.m
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -34,6 +34,7 @@ posky = ky>=0;
poskz = kz>=0;
kxeq0 = kx==0;
+kyeq0 = ky==0;
kzeq0 = kz==0;
omega = omega(posky,poskx,poskz);
@@ -56,6 +57,8 @@ subplot(1,nplots,iplot)
xlabel('$k_x$'); ylabel('$k_y$');
title('$\gamma(k_z=0)$');
end
+ if OPTIONS.INTERP; shading interp; end
+ caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzky
@@ -72,6 +75,8 @@ subplot(1,nplots,iplot)
xlabel('$k_z$'); ylabel('$k_y$');
title('$\gamma(k_x=0)$');
end
+ if OPTIONS.INTERP; shading interp; end
+ caxis(max(max(abs(toplot))).*[-1,1]);
iplot = iplot + 1;
end
if OPTIONS.kzkx
@@ -84,11 +89,12 @@ subplot(1,nplots,iplot)
title('$\max(\gamma)_{ky}$');
else
toplot = squeeze(real(omega(kyeq0,:,:)));
- pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
+ pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_z$'); ylabel('$k_x$');
title('$\gamma(k_y=0)$');
end
end
-shading interp
+if OPTIONS.INTERP; shading interp; end
+caxis(max(max(abs(toplot))).*[-1,1]);
colormap(bluewhitered);
end
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 03e36876..108efd96 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -17,41 +17,52 @@ end
FRAMES = unique(FRAMES);
t = DATA.Ts3D(FRAMES);
+
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
FIGURE.fig = figure; set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
for i = 1:2
- MODES_SELECTOR = i; %(2:Zonal, 1: NZonal, 3: ky=kx)
+ MODES_SELECTOR = i; %(2:Zonal, 1: NZonal)
if MODES_SELECTOR == 2
- if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES)))./max(abs(squeeze(x(1,ik,FRAMES)))),Nma);
- else
- plt = @(x,ik) movmean(abs(squeeze(x(1,ik,FRAMES))),Nma);
+ switch OPTIONS.ik
+ case 'sum'
+ plt_ = @(x,ik) movmean(squeeze(sum(abs((x(:,ik,FRAMES))),1)),Nma);
+ MODESTR = '$\sum_{k_y}$';
+ case 'max'
+ plt_ = @(x,ik) movmean(squeeze(max(abs((x(:,ik,FRAMES))),[],1)),Nma);
+ MODESTR = '$\max_{k_y}$';
+ otherwise
+ plt_ = @(x,ik) movmean(abs(squeeze(x(OPTIONS.ik,ik,FRAMES))),Nma);
+ MODESTR = ['$k_y=$',num2str(DATA.ky(OPTIONS.ik))];
end
kstr = 'k_x';
+ % Number max of modes to plot is kx>0 (1/2) of the non filtered modes (2/3)
+ Nmax = ceil(DATA.Nkx*1/3);
k = DATA.kx;
- MODESTR = 'Zonal modes';
elseif MODES_SELECTOR == 1
- if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES)))./max(abs(squeeze(x(ik,1,FRAMES)))),Nma);
- else
- plt = @(x,ik) movmean(abs(squeeze(x(ik,1,FRAMES))),Nma);
+ switch OPTIONS.ik
+ case 'sum'
+ plt_ = @(x,ik) movmean(squeeze(sum(abs((x(ik,:,FRAMES))),2)),Nma);
+ MODESTR = '$\sum_{k_x}$';
+ case 'max'
+ plt_ = @(x,ik) movmean(squeeze(max(abs((x(ik,:,FRAMES))),[],2)),Nma);
+ MODESTR = '$\max_{k_x}$';
+ otherwise
+ plt_ = @(x,ik) movmean(abs(squeeze(x(ik,OPTIONS.ik,FRAMES))),Nma);
+ MODESTR = ['$k_x=$',num2str(DATA.kx(OPTIONS.ik))];
end
kstr = 'k_y';
+ % Number max of modes to plot is ky>0 (1/1) of the non filtered modes (2/3)
+ Nmax = ceil(DATA.Nky*2/3);
k = DATA.ky;
- MODESTR = 'NZ modes';
- elseif MODES_SELECTOR == 3
- if NORMALIZED
- plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES)))./max(abs(squeeze(x(ik,ik,FRAMES)))),Nma);
- else
- plt = @(x,ik) movmean(abs(squeeze(x(ik,ik,FRAMES))),Nma);
- end
- kstr = 'k_y=k_x';
- k = DATA.ky;
- MODESTR = 'Diag modes';
end
+ if NORMALIZED
+ plt = @(x,ik) plt_(x,ik)./max(plt_(x,ik));
+ else
+ plt = @(x,ik) plt_(x,ik);
+ end
MODES = 1:numel(k);
% MODES = zeros(size(OPTIONS.K2PLOT));
@@ -71,7 +82,7 @@ for i = 1:2
amp(i_) = gr(2);
i_=i_+1;
end
- mod2plot = [2:OPTIONS.NMODES+1];
+ mod2plot = [2:min(Nmax,OPTIONS.NMODES+1)];
clr_ = jet(numel(mod2plot));
%plot
subplot(2,3,1+3*(i-1))
diff --git a/matlab/dbg_zBC_map.m b/matlab/dbg_zBC_map.m
index 48b1bffa..63de86cb 100644
--- a/matlab/dbg_zBC_map.m
+++ b/matlab/dbg_zBC_map.m
@@ -1,8 +1,19 @@
-Nkx = 16;
-Nky = 2;
+Nkx = 6;
+Nky = 3;
+
my = 0:(Nky-1);
mx = zeros(1,Nkx);
+PERIODIC = 1;
+Npol = 1;
+Nexc = 0;
+
+shear = 0.8;
+Ly = 120;
+Lx = 120;
+dky = 2*pi/Ly;
+
+
for ix = 1:Nkx
if(mod(Nkx,2) == 0)%even
mx_max = (Nkx/2);
@@ -22,12 +33,13 @@ for ix = 1:Nkx
end
disp(mx)
-Npol = 1;
-Nexc = 1;
-shear = 0.8;
-Ly = 120;
-dky = 2*pi/Ly;
-dkx = 2*pi*shear*dky/Nexc;
+
+if Nexc == 0 %% Adapt Nexc
+ dkx = 2*pi/Lx;
+ Nexc = ceil(0.9*2*pi*shear*dky/dkx);
+else
+ dkx = 2*pi*shear*dky/Nexc;
+end
kx = mx*dkx;
ky = my*dky;
@@ -38,12 +50,13 @@ for iy = 1:Nky
shift = 2*pi*shear*ky(iy)*Npol;
for ix = 1:Nkx
kx_shift = kx(ix) + shift;
- if 0%(kx_shift > kx_max)
+ if ((kx_shift > kx_max) && (~PERIODIC))
ikx_zBC_R(iy,ix) = nan;
else
ikx_zBC_R(iy,ix) = ix+(iy-1)*Nexc;
if(ikx_zBC_R(iy,ix) > Nkx)
- ikx_zBC_R(iy,ix) = ikx_zBC_R(iy,ix) - Nkx;
+% ikx_zBC_R(iy,ix) = ikx_zBC_R(iy,ix) - Nkx;
+ ikx_zBC_R(iy,ix) = mod(ikx_zBC_R(iy,ix)-1,Nkx)+1;
end
end
end
@@ -56,12 +69,13 @@ for iy = 1:Nky
shift = 2*pi*shear*ky(iy)*Npol;
for ix = 1:Nkx
kx_shift = kx(ix) - shift;
- if(kx_shift < kx_min)
+ if ((kx_shift < kx_min) && (~PERIODIC))
ikx_zBC_L(iy,ix) = nan;
else
ikx_zBC_L(iy,ix) = ix-(iy-1)*Nexc;
if(ikx_zBC_L(iy,ix) < 1)
- ikx_zBC_L(iy,ix) = ikx_zBC_L(iy,ix) + Nkx;
+% ikx_zBC_L(iy,ix) = ikx_zBC_L(iy,ix) + Nkx;
+ ikx_zBC_L(iy,ix) = mod(ikx_zBC_L(iy,ix)-1,Nkx)+1;
end
end
end
diff --git a/matlab/evolve_tracers.m b/matlab/evolve_tracers.m
index 4725ede3..5b7a4b40 100644
--- a/matlab/evolve_tracers.m
+++ b/matlab/evolve_tracers.m
@@ -2,18 +2,18 @@
SHOW_FILM = 1;
field2plot ='phi';
INIT = 'lin'; % lin (for a line)/ round (for a small round)/ gauss for random
-U_TIME = 1000; % >0 for frozen velocity at a given time, -1 for evolving field
+U_TIME = 500; % >0 for frozen velocity at a given time, -1 for evolving field
Evolve_U = 1; % 0 for frozen velocity at a given time, 1 for evolving field
-Tfin = 500;
+Tfin = 2000;
dt_ = 0.1;
Nstep = ceil(Tfin/dt_);
% Init tracers
-Np = 200; %number of tracers
+Np = 3; %number of tracers
% color = tcolors;
color = jet(Np);
tcolors = distinguishable_colors(Np); %Their colors
dimmed = 0; % to dimm the colormap in the background (infty = white, 0 normal color)
-Na = 10/dt_; %length of trace
+Na = 1000/dt_; %length of trace
Traj_x = zeros(Np,Nstep);
Traj_y = zeros(Np,Nstep);
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 3680bd89..705b7f3a 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
-tw = [00 1.3];
+tw = [00 40000];
fig = gcf;
axObjs = fig.Children;
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index af55df0c..1ccdb54b 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -56,7 +56,7 @@ for jt = 1:numel(DATA.Ts3D)
%
tmp = h5read([folder,momfile],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-
+%
tmp = h5read([folder,phifile],['/field/phi/',sprintf('%10.10d',it-1)]);
DATA.PHI(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index c4687314..ca11c1fd 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -11,6 +11,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
+ disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
% disp(['Q_x=',sprintf('%2.2e',Qx_infty_avg),'+-',sprintf('%2.2e',Qx_infty_std)]);
f_avg_z = squeeze(mean(DATA.PHI(:,:,:,:),3));
[~,ikzf] = max(squeeze(mean(abs(f_avg_z(1,:,its3D:ite3D)),3)));
@@ -62,7 +63,7 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
%% Figure
mvm = @(x) movmean(x,OPTIONS.NMVA);
- FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [500, 1000, 1000, 600])
subplot(311)
% yyaxis left
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index 6591df16..c0669f23 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -54,7 +54,7 @@ switch options.COMPXY
ynamecx = ['$\max_{k_x}',options.NAME,'$'];
ynamecy = ['$\max_{k_y}',options.NAME,'$'];
case 'zero'
- compx = @(x) x(:,data.Nx/2);
+ compx = @(x) x(:,1);
compy = @(x) x(1,:);
ynamecx = ['$',options.NAME,'(k_x=',num2str(data.kx(1)),')$'];
ynamecy = ['$',options.NAME,'(k_y=',num2str(data.ky(1)),')$'];
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 0804e892..baf1ab80 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -53,7 +53,7 @@ Lmin_ = min([Xmax_,Ymax_]);
Rx = Xmax_/Lmin_ * SCALE + (1-SCALE)*1.2;
Ry = Ymax_/Lmin_ * SCALE + (1-SCALE)*1.2;
ASPECT = [Rx, Ry, 1];
-DIMENSIONS = [100, 100, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
+DIMENSIONS = [500, 1000, OPTIONS.RESOLUTION*Rx, OPTIONS.RESOLUTION*Ry];
% Polar grid
POLARNAME = [];
if POLARPLOT
diff --git a/matlab/setup.m b/matlab/setup.m
index 71b95cdc..8e6245b7 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -21,6 +21,8 @@ GEOM.eps = EPS; % inverse aspect ratio
GEOM.kappa = KAPPA; % elongation
GEOM.delta = DELTA; % triangularity
GEOM.zeta = ZETA; % squareness
+GEOM.parallel_bc = ['''',PARALLEL_BC,''''];
+GEOM.shift_y = SHIFT_Y;
% Model parameters
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index adc09fc5..e708ca0f 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -33,6 +33,7 @@ fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
+fprintf(fid,[' parallel_bc = ', GEOM.parallel_bc,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
diff --git a/wk/analysis_gbms.m b/wk/analysis_gbms.m
index d1957270..376ae1b2 100644
--- a/wk/analysis_gbms.m
+++ b/wk/analysis_gbms.m
@@ -7,9 +7,9 @@
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/shearless_linear_cyclone/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/new_RH_test/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/RH_test_kine/';
-resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/KBM/';
+% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/KBM/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/TEM/';
-% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/ITG/';
+resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/ITG/';
% resdir = '/home/ahoffman/Documents/gbms/benchmark_HeLaZ/linear_cyclone/';
% resdir = '/home/ahoffman/molix/';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 48a6d160..5e609cbb 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -13,14 +13,17 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/LD_s_alpha_output_0.8/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
+folder = '/misc/gene_results/Z-pinch/HP_fig_2a_gyroLES/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_5e-2/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_gyroLES/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2c_mu_1e-2_muv_1e-1/';
% folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
-folder = '/misc/gene_results/Z-pinch/kN_2.0_HD_transport_spectrum_01/';
+% folder = '/misc/gene_results/Z-pinch/kN_2.0_HD_transport_spectrum_01/';
% folder = '/misc/gene_results/Z-pinch/kN_2.5_HD_transport_spectrum_01/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
@@ -33,6 +36,11 @@ folder = '/misc/gene_results/Z-pinch/kN_2.0_HD_transport_spectrum_01/';
% folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
% folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
% folder = '/misc/gene_results/miller/';
+% folder = '/misc/gene_results/CBC/KT_4.5_192x96x24x30x16_00/';
+% folder = '/misc/gene_results/CBC/KT_5.0_192x96x24x30x16_00/';
+% folder = '/misc/gene_results/CBC/KT_5.3_192x96x24x30x16_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
+
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -123,7 +131,7 @@ end
if 0
%% Show f_i(vpar,mu)
-options.times = 600:5000;
+options.times = 100:200;
options.specie = 'i';
options.PLT_FCT = 'contour';
options.folder = folder;
@@ -136,7 +144,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [600 5000];
+options.TIME = [4000 5000];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'n_i';
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index ad8a8c7b..f787c1c2 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -10,7 +10,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
%% Load the results
LOCALDIR = [gyacomodir,resdir,'/'];
-MISCDIR = ['/misc/',resdir,'/']; %For long term storage
+MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
@@ -34,7 +34,7 @@ disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.TAVG_0 = data.TJOB_SE(i_)+600;%0.4*data.Ts3D(end);
options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration
options.NCUT = 4; % Number of cuts for averaging and error estimation
-options.NMVA = 100; % Moving average for time traces
+options.NMVA = 50; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
options.INTERP = 0;
@@ -46,7 +46,7 @@ end
if 0
%% statistical transport averaging
gavg =[]; gstd = [];
-for i_ = 3:2:numel(data.TJOB_SE)
+for i_ = 1:2:numel(data.TJOB_SE)
% i_ = 3;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
@@ -75,7 +75,7 @@ options.PLAN = 'xy';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [000:1:8000];
+options.TIME = [6000:1:9000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -95,11 +95,12 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
-% options.NAME 'f_i';
-% options.PLAN = 'sx';
+options.PLAN = 'xy';
+options.NAME = 'f_i';
+options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [1000 3000 5000 8000 10000];
+options.TIME = [5 20 50 100 150];
+options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
@@ -127,11 +128,11 @@ if 0
options.SPAR = gene_data.vp';
options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [1500 5000];
+options.T = [100 150];
options.PLT_FCT = 'contour';
-options.ONED = 1;
+options.ONED = 0;
options.non_adiab = 0;
-options.SPECIE = 'e';
+options.SPECIE = 'i';
options.RMS = 1; % Root mean square i.e. sqrt(sum_k|f_k|^2) as in Gene
fig = plot_fa(data,options);
% save_figure(data,fig,'.png')
@@ -155,7 +156,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [2000 3000];
+options.TIME = [5000 9000];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
@@ -163,7 +164,7 @@ options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
-options.COMPXY = 'avg'; % avg/sum/max/zero/ 2D plot otherwise
+options.COMPXY = '2D'; % avg/sum/max/zero/ 2D plot otherwise
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
@@ -189,12 +190,14 @@ end
if 0
%% Mode evolution
-options.NORMALIZED = 0;
-options.K2PLOT = [0.1 0.5 1.0 2.0 6.0];
-options.TIME = [00:1200];
+options.NORMALIZED = 1;
+options.TIME = [6000:9000];
+options.KX_TW = [6000:9000]; %kx Z modes time window
+options.KY_TW = [6000:9000]; %ky Growth rate time window
options.NMA = 1;
-options.NMODES = 30;
-options.iz = 'avg';
+options.NMODES = 800;
+options.iz = 'avg'; % avg or index
+options.ik = 1; % sum, max or index
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
@@ -215,19 +218,25 @@ fig = plot_metric(data,options);
end
if 0
-%% linear growth rate for 3D fluxtube
-trange = [0 100];
-nplots = 1;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 2; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+options.GOK = 0; %plot 0: gamma 1: gamma/k 2: gamma^2/k^3
+lg = compute_fluxtube_growth_rate(data,options);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
end
if 0
%% linear growth rate for 3D Zpinch
-trange = [5 20];
+trange = [100 200];
options.keq0 = 1; % chose to plot planes at k=0 or max
options.kxky = 1;
-options.kzkx = 0;
+options.kzkx = 1;
options.kzky = 1;
+options.INTERP = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
% save_figure(data,fig,'.png')
end
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index 11cfa44a..6098397a 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -178,7 +178,7 @@ resdir ='';
% resdir ='Zpinch_rerun/Kn_2.0_256x64x9x5_Lx_240_Ly_120';
% resdir ='Zpinch_rerun/Kn_1.6_256x128x7x4';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6';
-resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
+% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
% resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_mu_0.5';
% resdir ='Zpinch_rerun/Kn_1.6_256x128x21x11';
@@ -234,9 +234,12 @@ resdir ='Zpinch_rerun/Kn_1.6_200x64x11x6_conv';
% resdir = 'Zpinch_rerun/DGGK_kN_1.6_200x64x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x5x3_nu_0.01';
-resdir = 'Zpinch_rerun/LRGK_kN_2.4_300x98x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/DGGK_kN_1.7_200x64x9x5_nu_0.01';
+% resdir = 'Zpinch_rerun/LRGK_kN_2.4_300x98x5x3_nu_0.01';
% resdir = 'Zpinch_rerun/LDGK_kN_1.9_200x64x5x3_nu_0.01';
+resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x5x3';
+% resdir = 'Zpinch_rerun/COLLESS_kN_1.7_200x64x9x5';
-JOBNUMMIN = 00; JOBNUMMAX = 03;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
resdir = ['results/',resdir];
run analysis_gyacomo
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 30e8b789..4dc2a135 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -31,11 +31,24 @@ gyacomodir = '/home/ahoffman/gyacomo/';
%% CBC
% resdir = 'CBC/64x32x16x5x3';
% resdir = 'CBC/64x128x16x5x3';
-% resdir = 'CBC/128x64x16x5x3';
+% resdir = 'CBC/old/128x64x16x5x3';
% resdir = 'CBC/96x96x16x3x2_Nexc_6';
-% resdir = 'CBC/128x96x16x3x2';
-resdir = 'CBC/192x96x16x3x2';
-% resdir = 'CBC/192x96x24x13x7';
+% resdir = 'CBC/128x96x16x3x2_Nexc_0';
+% resdir = 'CBC/old/192x96x24x13x7';
+
+% resdir = 'CBC/128x96x16x3x2_Nexc_0_periodic_chi';
+% resdir = 'CBC/64x32x16x3x2_Nexc_0_periodic_chi';
+% resdir = 'CBC/64x32x16x3x2_Nexc_0_Miller';
+% resdir = 'CBC/64x32x16x3x2_Nexc_1';
+% resdir = 'CBC/64x32x16x3x2_Nexc_2';
+% resdir = 'CBC/64x32x16x3x2_Nexc_3';
+% resdir = 'CBC/64x32x16x3x2_Nexc_3_change_dt';
+% resdir = 'CBC/64x32x16x3x2_Nexc_6';
+% resdir = 'CBC/64x32x16x3x2_Nexc_0';
+resdir = 'CBC/64x32x16x3x2_Nexc_0_change_dt';
+% resdir = 'CBC/128x96x16x3x2_Nexc_0';
+% resdir = 'CBC/128x96x16x3x2_Nexc_0';
+
% resdir = 'CBC/kT_11_128x64x16x5x3';
% resdir = 'CBC/kT_9_256x128x16x3x2';
% resdir = 'CBC/kT_4.5_128x64x16x13x3';
@@ -59,14 +72,14 @@ resdir = 'CBC/192x96x16x3x2';
% resdir = 'NL_KBM/192x64x24x5x3';
%% Linear CBC
% resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
+% resdir = 'lin_CBC/128x64x16_5x3_Lx_7.854_Ly_125.6637_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_6.96_nu_5.0e-02_DGDK_adiabe';
% resdir = 'testcases/miller_example';
% resdir = 'Miller/128x256x3x2_CBC_dt_5e-3';
%% CBC Miller
% resdir = 'GCM_CBC/daint/Miller_GX_comparison';
% resdir = 'GCM_CBC/daint/Salpha_GX_comparison';
%% RK3
-% resdir = 'dbg/SSPx_RK3_test';
-% resdir = 'dbg/SSPx_RK3_test/RK4';
-resdir = ['gyacomo_outputs/results/',resdir];
+% resdir = '';
+resdir = ['results/',resdir];
JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index e64bb2a7..e37ea07f 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -3,8 +3,8 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
-SIMID = 'lin_ITG'; % Name of the simulation
-% SIMID = 'dbg'; % Name of the simulation
+% SIMID = 'lin_ITG'; % Name of the simulation
+SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
@@ -34,23 +34,25 @@ PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 8; % real space x-gridpoints
-NY = 16; % '' y-gridpoints
+NX = 6; % real space x-gridpoints
+NY = 10; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.05; % Size of the squared frequency domain
+LY = 2*pi/0.1; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
%% GEOMETRY
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
-NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-EPS = 0.18; % inverse aspect ratio
+PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+SHIFT_Y = 0.0;
%% TIME PARMETERS
TMAX = 60; % Maximal time unit
DT = 1e-2; % Time step
@@ -103,11 +105,11 @@ setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
end
%% Load results
@@ -133,8 +135,8 @@ end
if 1
%% Ballooning plot
-options.time_2_plot = [120];
-options.kymodes = [0.3];
+options.time_2_plot = [10 50];
+options.kymodes = [0.1 0.2 0.4];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
--
GitLab