From a98b465e36f780a2e2c633c9e1a728d4c51a404f Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 25 Apr 2022 13:54:04 +0200
Subject: [PATCH] matlab scripts reorganization
---
matlab/compile_results.m | 2 +-
matlab/compute/LinearFit_s.m | 5 +-
.../compute/compute_3D_zpinch_growth_rate.m | 92 ++++++++++
{wk => matlab/compute}/mode_growth_meter.m | 0
{wk => matlab}/evolve_tracers.m | 0
{wk => matlab/load}/load_f_gene.m | 0
{wk => matlab/load}/quick_gene_load.m | 0
matlab/plot/real_plot_1D.m | 74 ++++++++-
matlab/setup.m | 10 +-
matlab/write_fort90.m | 12 +-
{wk => testcases}/Hallenbert.m | 0
.../benchmark_HeLaZ_gene_transport_zpinch.m | 0
{wk => testcases}/benchmark_HeLaZ_molix.m | 0
.../benchmark_linear_1D_entropy_mode.m | 0
{wk => testcases}/cyclone_test_case.m | 0
.../kobayashi_2015_linear_results.m | 0
{wk => testcases}/linear_1D_entropy_mode.m | 14 +-
{wk => testcases}/linear_damping.m | 0
wk/Dimits_shift_results.m | 95 -----------
wk/analysis_3D.m | 47 ++++--
wk/analysis_header.m | 16 +-
wk/gene_analysis_3D.m | 16 +-
wk/shearless_linear_fluxtube.m | 157 ------------------
23 files changed, 235 insertions(+), 305 deletions(-)
create mode 100644 matlab/compute/compute_3D_zpinch_growth_rate.m
rename {wk => matlab/compute}/mode_growth_meter.m (100%)
rename {wk => matlab}/evolve_tracers.m (100%)
rename {wk => matlab/load}/load_f_gene.m (100%)
rename {wk => matlab/load}/quick_gene_load.m (100%)
rename {wk => testcases}/Hallenbert.m (100%)
rename {wk => testcases}/benchmark_HeLaZ_gene_transport_zpinch.m (100%)
rename {wk => testcases}/benchmark_HeLaZ_molix.m (100%)
rename {wk => testcases}/benchmark_linear_1D_entropy_mode.m (100%)
rename {wk => testcases}/cyclone_test_case.m (100%)
rename {wk => testcases}/kobayashi_2015_linear_results.m (100%)
rename {wk => testcases}/linear_1D_entropy_mode.m (96%)
rename {wk => testcases}/linear_damping.m (100%)
delete mode 100644 wk/Dimits_shift_results.m
delete mode 100644 wk/shearless_linear_fluxtube.m
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 03aa07db..b98ac6c1 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -238,7 +238,7 @@ Nz = numel(z);
KPERP2 = KY.^2+KX.^2;
%% Add everything in output structure
% scaling
-DATA.scale = -(1/Nx/Ny)^2;
+DATA.scale = (1/Nx/Ny)^2;
% Fields
DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_;
DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
diff --git a/matlab/compute/LinearFit_s.m b/matlab/compute/LinearFit_s.m
index af790195..a2816ff8 100644
--- a/matlab/compute/LinearFit_s.m
+++ b/matlab/compute/LinearFit_s.m
@@ -30,5 +30,8 @@ function [gamma,t,gammaoft] = LinearFit_s(time,Na00abs)
% Return gamma(t) for amplitude ratio method
t = 0.5*(time(2:end)+time(1:end-1));
-
+
+ if isnan(gamma)
+ gamma =0;
+ end
end % ... end function
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
new file mode 100644
index 00000000..1cfeed82
--- /dev/null
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -0,0 +1,92 @@
+function [ omega, FIGURE ] = compute_3D_zpinch_growth_rate(DATA, TRANGE, OPTIONS)
+FIGURE.FIGNAME = ['growth_rate_kx=0_ky=0_planes',DATA.PARAMS];
+
+t = DATA.Ts3D;
+[~,its] = min(abs(t-TRANGE(1)));
+[~,ite] = min(abs(t-TRANGE(end)));
+nkx = DATA.Nkx; nky = DATA.Nky; nkz = DATA.Nz;
+% Remove AA part
+if DATA.Nx > 1
+ ikxnz = abs(DATA.PHI(:,1,1,1)) > 0;
+else
+ ikxnz = abs(DATA.PHI(:,1,1,1)) > -1;
+end
+ikynz = abs(DATA.PHI(1,:,1,1)) > 0;
+
+phi = fft(DATA.PHI(ikxnz,ikynz,:,:),[],3);
+
+omega = zeros(nkx,nky,nkz);
+
+for iz = 1:nkz
+ for iy = 1:nky
+ for ix = 1:nkx
+ omega(ix,iy,iz) = LinearFit_s(t(its:ite),squeeze(abs(phi(ix,iy,iz,its:ite))));
+ end
+ end
+end
+
+%% plot
+kx = DATA.kx; ky = DATA.ky; kz = [(0:nkz/2), (-nkz/2+1):-1];
+poskx = kx>=0;
+posky = ky>=0;
+poskz = kz>=0;
+
+kxeq0 = kx==0;
+kzeq0 = kz==0;
+
+omega = omega(poskx,posky,poskz);
+
+FIGURE.fig = figure;
+nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
+iplot = 1;
+
+if OPTIONS.kxky
+[Y_XY,X_XY] = meshgrid(ky(posky),kx(poskx));
+subplot(1,nplots,iplot)
+ if ~OPTIONS.keq0
+ toplot = squeeze(max(real(omega(:,:,:)),[],3));
+ pclr= pcolor(X_XY,Y_XY,toplot); set(pclr,'EdgeColor','none');
+ xlabel('$k_x$'); ylabel('$k_y$');
+ title('$\max(\gamma)_{kz}$');
+ else
+ toplot = squeeze(real(omega(:,:,kzeq0)));
+ pclr= pcolor(X_XY,Y_XY,toplot); set(pclr,'EdgeColor','none');
+ xlabel('$k_x$'); ylabel('$k_y$');
+ title('$\gamma(k_z=0)$');
+ end
+ iplot = iplot + 1;
+end
+if OPTIONS.kzky
+[Y_ZY,Z_ZY] = meshgrid(ky(posky),kz(poskz));
+subplot(1,nplots,iplot)
+ if ~OPTIONS.keq0
+ toplot = squeeze(max(real(omega(:,:,:)),[],1));
+ pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_x$'); ylabel('$k_y$');
+ title('$\max(\gamma)_{kx}$');
+ else
+ toplot = squeeze(real(omega(kxeq0,:,:)));
+ pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_z$'); ylabel('$k_y$');
+ title('$\gamma(k_x=0)$');
+ end
+ iplot = iplot + 1;
+end
+if OPTIONS.kzkx
+[X_ZX,Z_ZX] = meshgrid(kx(poskx),kz(poskz));
+subplot(1,nplots,iplot)
+ if ~OPTIONS.keq0
+ toplot = squeeze(max(real(omega(:,:,:)),[],2));
+ pclr= pcolor(Z_ZX,X_ZX,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_z$'); ylabel('$k_x$');
+ title('$\max(\gamma)_{ky}$');
+ else
+ toplot = squeeze(real(omega(:,kyeq0,:)));
+ pclr= pcolor(Z_ZY,Y_ZY,toplot'); set(pclr,'EdgeColor','none');
+ xlabel('$k_z$'); ylabel('$k_x$');
+ title('$\gamma(k_y=0)$');
+ end
+end
+shading interp
+colormap(bluewhitered);
+end
diff --git a/wk/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
similarity index 100%
rename from wk/mode_growth_meter.m
rename to matlab/compute/mode_growth_meter.m
diff --git a/wk/evolve_tracers.m b/matlab/evolve_tracers.m
similarity index 100%
rename from wk/evolve_tracers.m
rename to matlab/evolve_tracers.m
diff --git a/wk/load_f_gene.m b/matlab/load/load_f_gene.m
similarity index 100%
rename from wk/load_f_gene.m
rename to matlab/load/load_f_gene.m
diff --git a/wk/quick_gene_load.m b/matlab/load/quick_gene_load.m
similarity index 100%
rename from wk/quick_gene_load.m
rename to matlab/load/quick_gene_load.m
diff --git a/matlab/plot/real_plot_1D.m b/matlab/plot/real_plot_1D.m
index 5b7eda42..ec3217ff 100644
--- a/matlab/plot/real_plot_1D.m
+++ b/matlab/plot/real_plot_1D.m
@@ -31,22 +31,27 @@ switch options.NAME
yname = 'n_e';
end
-switch options.COMPXY
+switch options.COMPX
case 'avg'
compx = @(x) mean(x,1);
- compy = @(x) mean(x,2);
ynamex = ['$\langle ',yname,'\rangle_x$'];
- ynamey = ['$\langle ',yname,'\rangle_y$'];
otherwise
compx = @(x) x(1,:);
+ ynamex = ['$',yname,'(y,x=0)$'];
+end
+
+switch options.COMPY
+ case 'avg'
+ compy = @(x) mean(x,2);
+ ynamey = ['$\langle ',yname,'\rangle_y$'];
+ otherwise
compy = @(x) x(:,1);
- ynamex = ['$',yname,'(x=0)$'];
- ynamey = ['$',yname,'(y=0)$'];
+ ynamey = ['$',yname,'(x,y=0)$'];
end
-set(gcf, 'Position', [20 50 5000 2000])
-subplot(1,2,1)
+set(gcf, 'Position', [20 50 1000 500])
+subplot(1,3,1)
X = data.x;
xname = '$x$';
@@ -76,7 +81,7 @@ subplot(1,2,1)
xlabel(xname); ylabel(ynamey)
xlim([min(X),max(X)]);
-subplot(1,2,2)
+subplot(1,3,2)
X = data.y;
xname = '$y$';
@@ -108,6 +113,59 @@ subplot(1,2,2)
legend('show','Location','eastoutside');
xlabel(xname); ylabel(ynamex)
xlim([min(X),max(X)]);
+
+%% Z plot
+if data.Nz > 1
+ options.PLAN = 'yz';
+ options.COMP = options.COMPX;
+ options.POLARPLOT = 0;
+ options.INTERP = 0;
+
+ toplot = process_field(data,options);
+ t = data.Ts3D; frames = toplot.FRAMES;
+end
+
+switch options.COMPZ
+ case 'avg'
+ compz = @(x) mean(x,1);
+ ynamez = ['$\langle ',yname,'\rangle_y$'];
+ otherwise
+ compz = @(x) x(1,:);
+ ynamez = ['$',yname,'(z,y=0)$'];
+end
+
+subplot(1,3,3)
+
+ X = data.z;
+ xname = '$z$';
+ switch options.COMPT
+ case 'avg'
+ Y = compz(mean(toplot.FIELD(:,:,:),3));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t\in[',num2str(t(frames(1))),',',num2str(t(frames(end))),']$'],...
+ 'Color',colors(1,:)); hold on;
+ legend('show')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Y = compx(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = squeeze(Y);
+ if options.NORM
+ Y = Y./max(abs(Y));
+ end
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
+ legend('show','Location','eastoutside')
+ end
+
+ grid on
+% title('HeLaZ $k_y$ transport spectrum');
+ legend('show','Location','eastoutside');
+ xlabel(xname); ylabel(ynamez)
+ xlim([min(X),max(X)]);
end
diff --git a/matlab/setup.m b/matlab/setup.m
index d1afb57a..00269b32 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -10,10 +10,12 @@ GRID.Lx = LX; % x length
GRID.Ny = NY; % y ''
GRID.Ly = LY; % y ''
GRID.Nz = NZ; % z resolution
-GRID.q0 = Q0; % q factor
-GRID.shear = SHEAR; % shear
-GRID.eps = EPS; % inverse aspect ratio
if SG; GRID.SG = '.true.'; else; GRID.SG = '.false.';end;
+% Geometry
+GEOM.geom = ['''',GEOMETRY,''''];
+GEOM.q0 = Q0; % q factor
+GEOM.shear = SHEAR; % shear
+GEOM.eps = EPS; % inverse aspect ratio
% Model parameters
MODEL.CLOS = CLOS;
MODEL.NL_CLOS = NL_CLOS;
@@ -159,7 +161,7 @@ if ~exist(BASIC.MISCDIR, 'dir')
mkdir(BASIC.MISCDIR)
end
%% Compile and WRITE input file
-INPUT = write_fort90(OUTPUTS,GRID,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC);
+INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
MAKE = 'cd ..; make; cd wk';
% system(MAKE);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 574b6d02..d35641f0 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -1,4 +1,4 @@
-function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC)
+function [INPUT] = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC)
% Write the input script "fort.90" with desired parameters
INPUT = ['fort_',sprintf('%2.2d',OUTPUTS.job2load+1),'.90'];
fid = fopen(INPUT,'wt');
@@ -20,12 +20,16 @@ fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']);
fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']);
fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
-fprintf(fid,[' q0 = ', num2str(GRID.q0),'\n']);
-fprintf(fid,[' shear = ', num2str(GRID.shear),'\n']);
-fprintf(fid,[' eps = ', num2str(GRID.eps),'\n']);
fprintf(fid,[' SG = ', GRID.SG,'\n']);
fprintf(fid,'/\n');
+fprintf(fid,'&GEOMETRY\n');
+fprintf(fid,[' geom = ', GEOM.geom,'\n']);
+fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']);
+fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']);
+fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
+fprintf(fid,'/\n');
+
fprintf(fid,'&OUTPUT_PAR\n');
fprintf(fid,[' nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
fprintf(fid,[' nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
diff --git a/wk/Hallenbert.m b/testcases/Hallenbert.m
similarity index 100%
rename from wk/Hallenbert.m
rename to testcases/Hallenbert.m
diff --git a/wk/benchmark_HeLaZ_gene_transport_zpinch.m b/testcases/benchmark_HeLaZ_gene_transport_zpinch.m
similarity index 100%
rename from wk/benchmark_HeLaZ_gene_transport_zpinch.m
rename to testcases/benchmark_HeLaZ_gene_transport_zpinch.m
diff --git a/wk/benchmark_HeLaZ_molix.m b/testcases/benchmark_HeLaZ_molix.m
similarity index 100%
rename from wk/benchmark_HeLaZ_molix.m
rename to testcases/benchmark_HeLaZ_molix.m
diff --git a/wk/benchmark_linear_1D_entropy_mode.m b/testcases/benchmark_linear_1D_entropy_mode.m
similarity index 100%
rename from wk/benchmark_linear_1D_entropy_mode.m
rename to testcases/benchmark_linear_1D_entropy_mode.m
diff --git a/wk/cyclone_test_case.m b/testcases/cyclone_test_case.m
similarity index 100%
rename from wk/cyclone_test_case.m
rename to testcases/cyclone_test_case.m
diff --git a/wk/kobayashi_2015_linear_results.m b/testcases/kobayashi_2015_linear_results.m
similarity index 100%
rename from wk/kobayashi_2015_linear_results.m
rename to testcases/kobayashi_2015_linear_results.m
diff --git a/wk/linear_1D_entropy_mode.m b/testcases/linear_1D_entropy_mode.m
similarity index 96%
rename from wk/linear_1D_entropy_mode.m
rename to testcases/linear_1D_entropy_mode.m
index 1a74761d..f32cb24e 100644
--- a/wk/linear_1D_entropy_mode.m
+++ b/testcases/linear_1D_entropy_mode.m
@@ -6,9 +6,9 @@ default_plots_options
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 0.01; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_N = 1.6; % Density gradient drive
+K_N = 2.0; % 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)
@@ -18,12 +18,14 @@ NY = 1; % '' y-gridpoints
LX = 120; % Size of the squared frequency domain
LY = 1; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
+SG = 1; % Staggered z grids option
+%% GEOMETRY
+GEOMETRY= 'Z-pinch';
Q0 = 1.0; % safety factor
SHEAR = 0.0; % magnetic shear
EPS = 0.0; % inverse aspect ratio
-SG = 1; % Staggered z grids option
%% TIME PARMETERS
-TMAX = 200; % Maximal time unit
+TMAX = 100; % 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
@@ -37,7 +39,7 @@ 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 = 'LD';
+CO = 'DG';
GKCO = 0; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -99,7 +101,7 @@ for i = 1:Nparam
system(['rm fort*.90']);
% Run linear simulation
if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz3 1 1 1 0; cd ../../../wk'])
% disp([param_name,'=',num2str(K_N)]);
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz3 1 6 0 > out.txt; cd ../../../wk']);
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ./../../../bin/helaz 1 2 0; cd ../../../wk'])
diff --git a/wk/linear_damping.m b/testcases/linear_damping.m
similarity index 100%
rename from wk/linear_damping.m
rename to testcases/linear_damping.m
diff --git a/wk/Dimits_shift_results.m b/wk/Dimits_shift_results.m
deleted file mode 100644
index 3f846c77..00000000
--- a/wk/Dimits_shift_results.m
+++ /dev/null
@@ -1,95 +0,0 @@
-%% Nu = 1e-2
-KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.5];
-Ginf_LDGK = [5.5e-2 1.3e-1 2.7e-1 4.6e-1 1.0e+0 0.0e+0 0.0e+0];
-conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-Ginf_SGGK = [1.2e-3 2.8e-3 1.7e-2 8.9e-2 3.0e+0 0.0e+0 0.0e+0];
-conv_SGGK = [9.0e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-
-Ginf_DGGK = [1.1e-4 5.6e-4 9.0e-4 6.5e-3 7.0e+0 0.0e+0 0.0e+0];
-conv_DG = [1.6e-4 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-figure
-semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
-
-
-semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,conv_SGGK,'x-k','DisplayName','conv'); hold on;
-
-
-semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_DG,'x-k','DisplayName','conv'); hold on;
-
-title('$\nu=0.01$');
-xlabel('Drive');
-ylabel('Transport');
-
-%% Nu = 5e-2
-KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LDGK = [0.0e+0 0.0e+0 0.0e+0 8.3e-1 0.0e+0];
-conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-Ginf_SGGK = [0.0e+0 0.0e+0 0.0e+0 2.0e-1 0.0e+0];
-conv_SGGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-Ginf_SGDK = [0.0e+0 0.0e+0 0.0e+0 5.6e-3 0.0e+0];
-conv_SGDK = [0.0e+0 0.0e+0 0.0e+0 2.0e-3 0.0e+0];
-
-Ginf_DGGK = [0.0e+0 0.0e+0 0.0e+0 1.5e-2 0.0e+0];
-figure
-semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau GK'); hold on;
-semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama GK'); hold on;
-semilogy(KN,Ginf_SGDK,'o-c','DisplayName','Sugama DK'); hold on;
-semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty GK'); hold on;
-semilogy(KN,conv_LDGK,'x--g','DisplayName','changing params'); hold on;
-semilogy(KN,conv_SGGK,'x--b','DisplayName','changing params'); hold on;
-semilogy(KN,conv_SGDK,'x--c','DisplayName','changing params'); hold on;
-title('$\nu=0.05$');
-xlabel('Drive');
-ylabel('Transport');
-
-
-%% Nu = 1e-1
-KN = [ 1.5 1.6 1.7 1.8 1.9 2.0 2.5];
-Ginf_LDGK = [2.5e-2 2.5e-1 6.3e-1 1.3e+0 2.3e+0 4.3e+0 0.0e+0];
-conv_LDGK = [0.0e+0 2.5e-1 0.0e+0 1.2e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_SGGK = [0.0e-9 1.3e-2 9.1e-2 3.0e-1 7.8e-1 1.4e+0 3.5e+1];
-Ginf_DGGK = [2.0e-4 6.4e-3 3.6e-2 3.1e-1 0.0e+0 3.0e-1 3.3e+1];
-
-figure
-semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_LDGK,'x--k','DisplayName','changing params'); hold on;
-title('$\nu=0.1$');
-xlabel('Drive');
-ylabel('Transport');
-
-
-%% Nu = 5e-1
-KN = [ 1.5 1.6 1.7 1.8 1.9];
-Ginf_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-conv_LDGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-Ginf_SGGK = [0.0e-9 3.9e-2 2.7e-1 6.6e-1 1.6e+0];
-conv_SGGK = [0.0e-9 0.0e+0 3.0e-1 0.0e+0 1.6e+0];
-Ginf_DGGK = [0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-
-figure
-semilogy(KN,Ginf_LDGK,'d--g','DisplayName','Landau'); hold on;
-semilogy(KN,Ginf_SGGK,'o-b','DisplayName','Sugama'); hold on;
-semilogy(KN,Ginf_DGGK,'^-r','DisplayName','Dougherty'); hold on;
-semilogy(KN,conv_LDGK,'x--k','DisplayName','changing params'); hold on;
-semilogy(KN,conv_SGGK,'x--k','DisplayName','changing params'); hold on;
-title('$\nu=0.5$');
-xlabel('Drive');
-ylabel('Transport');
-
-%% Primary mode stability region (confirmed with GENE)
-
-stable = [1.5 0.1; 1.5 0.05];
-unstable = [1.5 0.01];
-
-figure
-semilogy(stable(:,1),stable(:,2),'og'); hold on;
-semilogy(unstable(:,1),unstable(:,2),'or');
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index c7a82f47..9608760e 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,4 +1,3 @@
-helazdir = '/home/ahoffman/HeLaZ/';
addpath(genpath([helazdir,'matlab'])) % ... add
addpath(genpath([helazdir,'matlab/plot'])) % ... add
addpath(genpath([helazdir,'matlab/compute'])) % ... add
@@ -11,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 = 0;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -20,9 +19,9 @@ default_plots_options
disp('Plots')
FMT = '.fig';
-if 0
+if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 0; TAVG_1 = 400; % Averaging times duration
+TAVG_0 = 300; TAVG_1 = 400; % Averaging times duration
compz = 'avg';
% chose your field to plot in spacetime diag (uzf,szf,Gx)
fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi',1,compz);
@@ -46,12 +45,12 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'yz';
+options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
-options.COMP = 1;
+options.COMP = 9;%'avg';
% options.TIME = dat.Ts5D;
-options.TIME = 0:0.5:100;
+options.TIME = 800:820;
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -63,17 +62,17 @@ if 0
options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
-options.NAME = '\phi';
-% options.NAME = 'v_y';
+% options.NAME = '\phi';
+% options.NAME = 'n_i';
% options.NAME = 'N_i^{00}';
-% options.NAME = 'T_i';
+options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'yz';
+options.PLAN = 'xy';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [10 20 100];
+options.TIME = [300 350 400];
data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -134,14 +133,15 @@ end
if 0
%% 1D real plot
-options.TIME = 1500:100:2500;
+options.TIME = 20;
options.NORM = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME ='\Gamma_x';
-options.NAME ='s_y';
+% options.NAME ='s_y';
+options.COMPX = 1;
+options.COMPY = 1;
options.COMPZ = 1;
-options.COMPXY = 'avg';
options.COMPT = 'avg';
options.MOVMT = 1;
fig = real_plot_1D(data,options);
@@ -152,7 +152,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 1;
options.K2PLOT = 0.01:0.01:1.0;
-options.TIME = 20:1:600;
+options.TIME = 20:1:60;
options.NMA = 1;
options.NMODES = 20;
fig = mode_growth_meter(data,options);
@@ -172,9 +172,20 @@ if 0
fig = plot_metric(data);
end
-if 1
+if 0
%% linear growth rate for 3D fluxtube
trange = [0 100];
nplots = 1;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
+end
+
+if 0
+%% linear growth rate for 3D Zpinch
+trange = [30 40];
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 1;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
end
\ No newline at end of file
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 657b52a8..c9270219 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -1,8 +1,16 @@
-%% Directory of the simulation
+% Directory of the code
+helazdir = '/home/ahoffman/HeLaZ/';
+% Directory of the simulation
% if 1% Local results
outfile ='';
-outfile ='linear_shearless_cyclone/dbg';
+% outfile ='3D_Zpinch/dbg';
+% outfile ='3D_Zpinch/entropy_mode';
+% outfile ='3D_Zpinch/ITG_Npol_gt_1';
+% outfile ='3D_Zpinch/Ivanov_ITG';
+% outfile ='3D_Zpinch/NL_ITG';
% outfile ='debug/test_zpar/';
-outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
+% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_0.8/';
-analysis_3D
\ No newline at end of file
+outfile ='shearless_cyclone/marconi/';
+% outfile ='shearless_cyclone/180x180x20x4x2_L_120_CBC_0.8/';
+call analysis_3D
\ No newline at end of file
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
index 2289b64f..cfde12ed 100644
--- a/wk/gene_analysis_3D.m
+++ b/wk/gene_analysis_3D.m
@@ -1,7 +1,9 @@
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
-folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.8/';
+% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.8/';
+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/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
gene_data = load_gene_data(folder);
@@ -20,20 +22,20 @@ 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_i';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'xz';
% options.NAME ='f_e';
% options.PLAN = 'sx';
-options.COMP = 9;
-options.TIME = [100 200 400];
-data.a = data.EPS * 2000;
+options.COMP = 'avg';
+options.TIME = [100 300 900];
+gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig)
end
diff --git a/wk/shearless_linear_fluxtube.m b/wk/shearless_linear_fluxtube.m
deleted file mode 100644
index 51ed3e45..00000000
--- a/wk/shearless_linear_fluxtube.m
+++ /dev/null
@@ -1,157 +0,0 @@
-RUN = 1; % To run or just to load
-addpath(genpath('../matlab')) % ... add
-default_plots_options
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Set Up parameters
-CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% PHYSICAL PARAMETERS
-NU = 0.0; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-K_N = 2.22; % Density gradient drive
-K_T = 6.0; % Temperature '''
-SIGMA_E = 0.05196; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-KIN_E = 1; % Kinetic (1) or adiabatic (0) electron model
-%% GRID PARAMETERS
-NX = 2; % real space x-gridpoints
-NY = 16; % '' y-gridpoints
-LX = 2*pi/0.25; % Size of the squared frequency domain
-LY = 2*pi/0.25; % Size of the squared frequency domain
-NZ = 24; % number of perpendicular planes (parallel grid)
-Q0 = 2.7; % safety factor
-SHEAR = 0.0; % magnetic shear
-EPS = 0.18; % inverse aspect ratio
-SG = 1; % Staggered z grids option
-%% TIME PARMETERS
-TMAX = 40; % Maximal time unit
-DT = 8e-3; % Time step
-SPS0D = 1; % Sampling per time unit for 2D arrays
-SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 10; % Sampling per time unit for 2D arrays
-SPS5D = 1/100; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-JOB2LOAD= -1;
-%% OPTIONS
-SIMID = 'shearless_fluxtube'; % Name of the simulation
-% 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: gyrofluid closure (p+2j<=Pmax))
-NL_CLOS =-1; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
-KERN = 0; % Kernel model (0 : GK)
-%% OUTPUTS
-W_DOUBLE = 0;
-W_GAMMA = 1; W_HF = 1;
-W_PHI = 1; W_NA00 = 1;
-W_DENS = 1; W_TEMP = 1;
-W_NAPJ = 1; W_SAPJ = 0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% unused
-HD_CO = 0.5; % Hyper diffusivity cutoff ratio
-kmax = NX*pi/LX;% Highest fourier mode
-NU_HYP = 0.0; % Hyperdiffusivity coefficient
-MU = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
-MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
-K_E = 0.00; % Electrostat '''
-GRADB = 1.0;
-CURVB = 1.0;
-INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0; INIT_PHI= 0;
-NOISE0 = 0.0e-4; % Init noise amplitude
-BCKGD0 = 1.0e-4; % Init background
-LAMBDAD = 0.0;
-KXEQ0 = 0; % put kx = 0
-LINEARITY = 0; % activate non-linearity (is cancelled if KXEQ0 = 1)
-%% PARAMETER SCANS
-
-%% Parameter scan over PJ
-% PA = [2 4];
-% JA = [1 2];
-PA = [4];
-JA = [2];
-DTA= DT*ones(size(JA));%./sqrt(JA);
-% DTA= DT;
-mup_ = MU_P;
-muj_ = MU_J;
-Nparam = numel(PA);
-param_name = 'PJ';
-Nkx = NX;
-Nky = NY/2+1;
-gamma_Ni00 = zeros(Nparam,Nkx,Nky);
-gamma_Nipj = zeros(Nparam,Nkx,Nky);
-gamma_phi = zeros(Nparam,Nkx,Nky);
-for i = 1:Nparam
- % Change scan parameter
- PMAXE = PA(i); PMAXI = PA(i);
- JMAXE = JA(i); JMAXI = JA(i);
- DT = DTA(i);
- setup
-% system(['rm fort*.90']);
- % Run linear simulation
- if 0
- system(['cd ../results/',SIMID,'/',PARAMS,'/; ./../../../bin/helaz 0; cd ../../../wk'])
- end
-% Load and process results
- %%
- filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
- load_results
-
- for ikx = 1:Nkx
- for iky = 1:Nky
- tend = max(Ts3D(abs(Ni00(ikx,iky,1,:))~=0));
- tstart = 0.2*tend;
- [~,itstart] = min(abs(Ts3D-tstart));
- [~,itend] = min(abs(Ts3D-tend));
- trange = itstart:itend;
- % exp fit on moment 00
- X_ = Ts3D(trange); Y_ = squeeze(abs(Ni00(ikx,iky,1,trange)));
- gamma_Ni00(i,ikx,iky) = LinearFit_s(X_,Y_);
- % exp fit on phi
- X_ = Ts3D(trange); Y_ = squeeze(abs(PHI(ikx,iky,1,trange)));
- gamma_phi (i,ikx,iky) = LinearFit_s(X_,Y_);
- end
- 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));
- if 0
- %% Fit verification
- figure;
- for i = 1:1:Nky
- X_ = Ts3D(:); Y_ = squeeze(abs(Ni00(1,i,1,:)));
- semilogy(X_,Y_,'DisplayName',['k=',num2str(ky(i))]); hold on;
- end
- end
-end
-
-if 1
-%% Plot
-SCALE = 1;%sqrt(2);
-fig = figure; FIGNAME = 'linear_study';
-plt = @(x) squeeze(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);
- plot(plt(SCALE*ky(1:Nky)),plt(gamma_phi(i,2,1:end)),...
- 'Color',clr,...
- 'LineStyle',linestyle{1},'Marker','^',...
-...% 'DisplayName',['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, $P=',num2str(PA(i)),'$, $J=',num2str(JA(i)),'$']);
- 'DisplayName',[CONAME,', $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$']);
- hold on;
- end
- grid on; xlabel('$k_y\rho_s^{R}$'); ylabel('$\gamma(\phi)L_\perp/c_s$'); xlim([0.0,max(ky)]);
- title(['$\kappa_N=',num2str(K_N),'$, $\nu_{',CONAME,'}=',num2str(NU),'$'])
- legend('show'); %xlim([0.01,10])
-saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.fig']);
-saveas(fig,[SIMDIR,'/',PARAMS,'/gamma_vs_',param_name,'_',PARAMS,'.png']);
-end
-
-if 0
-%% Trajectories of some modes
-figure;
-for i = 1:1:Nky
- semilogy(Ts3D,squeeze(abs(Ni00(1,i,1,:))),'DisplayName',['k=',num2str(ky(i))]); hold on;
-end
-end
\ No newline at end of file
--
GitLab