From 7e6cccbcc49d49e6d4d26cd631b4839540a397be Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 4 Nov 2022 15:45:52 +0100
Subject: [PATCH] Scripts save
---
.../compute/compute_3D_zpinch_growth_rate.m | 12 +-
matlab/extract_fig_data.m | 4 +-
matlab/plot/real_plot_1D.m | 2 +-
matlab/plot/statistical_transport_averaging.m | 2 +-
matlab/setup.m | 4 +-
scripts/scan_kN_pipeline.sh | 59 ++++++
testcases/zpinch_example/fort.90 | 18 +-
wk/CBC_nu_CO_scan.m | 172 ++++++++++------
wk/Zpinch_coll_scan_kN_1.7.m | 127 +++++++++++-
wk/analysis_gene.m | 37 ++--
wk/analysis_gyacomo.m | 51 ++---
wk/header_2DZP_results.m | 24 ++-
wk/lin_3D_Zpinch.m | 194 ++++++++++++++++++
wk/{lin_EPY.m => lin_ETPY.m} | 32 +--
wk/lin_ITG.m | 13 +-
15 files changed, 598 insertions(+), 153 deletions(-)
create mode 100644 scripts/scan_kN_pipeline.sh
create mode 100644 wk/lin_3D_Zpinch.m
rename wk/{lin_EPY.m => lin_ETPY.m} (87%)
diff --git a/matlab/compute/compute_3D_zpinch_growth_rate.m b/matlab/compute/compute_3D_zpinch_growth_rate.m
index 3a053fa6..ce129abd 100644
--- a/matlab/compute/compute_3D_zpinch_growth_rate.m
+++ b/matlab/compute/compute_3D_zpinch_growth_rate.m
@@ -13,6 +13,8 @@ else
end
ikynz = abs(DATA.PHI(:,1,1,1)) > 0;
+ikynz(1) = 1; ikxnz(1) = 1; %put k=0 in the analysis
+
phi = fft(DATA.PHI(ikynz,ikxnz,:,:),[],3);
omega = zeros(nky,nkx,nkz);
@@ -41,16 +43,16 @@ nplots = OPTIONS.kxky + OPTIONS.kzkx + OPTIONS.kzky;
iplot = 1;
if OPTIONS.kxky
-[X_XY,Y_XY] = meshgrid(ky(posky),kx(poskx));
+[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');
+ 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');
+ pclr= pcolor(X_XY,Y_XY,toplot'); set(pclr,'EdgeColor','none');
xlabel('$k_x$'); ylabel('$k_y$');
title('$\gamma(k_z=0)$');
end
@@ -65,7 +67,7 @@ subplot(1,nplots,iplot)
xlabel('$k_x$'); ylabel('$k_y$');
title('$\max(\gamma)_{kx}$');
else
- toplot = squeeze(real(omega(kxeq0,:,:)));
+ 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)$');
@@ -81,7 +83,7 @@ subplot(1,nplots,iplot)
xlabel('$k_z$'); ylabel('$k_x$');
title('$\max(\gamma)_{ky}$');
else
- toplot = squeeze(real(omega(:,kyeq0,:)));
+ 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)$');
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index a036d184..3680bd89 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 5000];
+tw = [00 1.3];
fig = gcf;
axObjs = fig.Children;
@@ -22,7 +22,7 @@ end
mvm = @(x) movmean(x,1);
shift = 0;%X_(n0);
% shift = 0;
-skip = 50;
+skip = 1;
figure;
plot(mvm(X_(n0:skip:n1)-shift),mvm(Y_(n0:skip:n1))); hold on;
diff --git a/matlab/plot/real_plot_1D.m b/matlab/plot/real_plot_1D.m
index d0e93cd5..8bdc5464 100644
--- a/matlab/plot/real_plot_1D.m
+++ b/matlab/plot/real_plot_1D.m
@@ -154,7 +154,7 @@ subplot(1,nplots,3)
legend('show')
otherwise
for it = 1:numel(toplot.FRAMES)
- Y = compx(toplot.FIELD(:,:,toplot.FRAMES(it)));
+ Y = compx(toplot.FIELD(:,:,it));
Y = squeeze(Y);
if options.NORM
Y = Y./max(abs(Y));
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index 3a7cc924..ffc0395e 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -1,5 +1,5 @@
-function [ fig ] = statistical_transport_averaging( data, options )
+function [ fig, Gx_infty_avg, Gx_infty_std ] = statistical_transport_averaging( data, options )
scale = data.scale;
Trange = options.T;
[~,it0] = min(abs(Trange(1)-data.Ts0D));
diff --git a/matlab/setup.m b/matlab/setup.m
index 3b8be6c7..cf9e1f15 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -66,10 +66,10 @@ switch CO
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_50_kpm_4.0.h5''';
+ COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
- COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
+% COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
end
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
diff --git a/scripts/scan_kN_pipeline.sh b/scripts/scan_kN_pipeline.sh
new file mode 100644
index 00000000..8b9cef7e
--- /dev/null
+++ b/scripts/scan_kN_pipeline.sh
@@ -0,0 +1,59 @@
+#! /bin/bash
+# lastjid=$(sbatch submit_00.cmd)
+
+nu_=0.01
+dnu=0.0
+kn_=1.7
+dkn=0.2
+kt_=0.425
+dkt=0.05
+Lx_=150
+dLx=030
+
+Tm_=4000
+dTm=1000
+
+# First submit
+istart=0
+lastji=$(sbatch submit_0$istart.cmd)
+lastjid=${lastji##* }
+echo $lastji
+
+for i in {1..5}; do
+ # Setup indices of job id (current and previous one)
+ im1=$(awk "BEGIN {print $i-1}")
+ idm1=$(printf '%02d' $im1)
+ id=$(printf '%02d' $i)
+
+ # Create new submit file from older one
+ awk -v "ID=$id" '{
+ if (NR == 8) print "#SBATCH --error=err_"ID".txt";
+ else if (NR == 9) print "#SBATCH --output=out_"ID".txt";
+ else if (NR == 12) print "srun --cpu-bind=cores ./gyacomo 2 24 1 "ID;
+ else print $0}' submit_$idm1.cmd > submit_$id.cmd
+
+ # Create new fort file from older one
+ awk -v "NU=$nu_" -v "TM=$Tm_" -v "J2L=$im1" -v "KN=$kn_" -v "KT=$kt_" -v "LX=$Lx_" '{
+ if (NR == 04) print " tmax = "TM;
+ else if (NR == 40) print " job2load = "J2L;
+ else if (NR == 13) print " Lx = "LX;
+ else if (NR == 54) print " nu = "NU;
+ else if (NR == 61) print " K_Ne = "KN;
+ else if (NR == 62) print " K_Ni = "KN;
+ else if (NR == 63) print " K_Te = "KT;
+ else if (NR == 64) print " K_Ti = "KT;
+ else print $0}' fort_$idm1.90 > fort_$id.90
+
+ # Retrieve last jobid and launch next job with dep
+ lastji=$(sbatch --dependency=afterok:$lastjid submit_$id.cmd)
+ lastjid=${lastji##* }
+ echo $lastjid
+
+ # Increment variables
+ Lx_=$(awk "BEGIN {print $Lx_+$dLx}")
+ nu_=$(awk "BEGIN {print $nu_+$dnu}")
+ kn_=$(awk "BEGIN {print $kn_+$dkn}")
+ kt_=$(awk "BEGIN {print $kt_+$dkt}")
+ Tm_=$(awk "BEGIN {print $Tm_+$dTm}")
+done
+
diff --git a/testcases/zpinch_example/fort.90 b/testcases/zpinch_example/fort.90
index f769f6b8..53a1a749 100644
--- a/testcases/zpinch_example/fort.90
+++ b/testcases/zpinch_example/fort.90
@@ -1,7 +1,7 @@
&BASIC
nrun = 100000000
- dt = 0.005
- tmax = 100
+ dt = 0.01
+ tmax = 200
maxruntime = 356400
/
&GRID
@@ -58,20 +58,20 @@
sigma_i = 1
q_e = -1
q_i = 1
- K_Ne = 2.22
- K_Te = 6.96
- K_Ni = 2.22
- K_Ti = 6.96
+ K_Ne = 2.0
+ K_Te = 0.4
+ K_Ni = 2.0
+ K_Ti = 0.4
GradB = 1
CurvB = 1
lambdaD = 0
- beta = 1e-4
+ beta = 0
/
&COLLISION_PAR
- collision_model = 'SG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
gyrokin_CO = .false.
interspecies = .true.
- mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+ !mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
INIT_OPT = 'phi'
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
index d14d69ba..7f69c607 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_CO_scan.m
@@ -1,47 +1,61 @@
-% NU_a = [0.05 0.15 0.25 0.35 0.45];
-NU_a = [0:0.1:0.5];
-g_max= NU_a*0;
-g_avg= NU_a*0;
-g_std= NU_a*0;
-k_max= NU_a*0;
-CO = 'DG';
-
-K_T = 7;
-DT = 5e-3;
-TMAX = 20;
-ky_ = 0.3;
-SIMID = 'linear_CBC_nu_scan_kT_7_ky_0.3_DGGK'; % Name of the simulation
+gyacomodir = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_dbg';
+EXECNAME = 'gyacomo';
+%%
+NU_a = [0.005 0.01:0.01:0.1];
+% P_a = [2 4 6 8 10 12 16];
+% NU_a = 0.05;
+P_a = [10];
+
+
+CO = 'SG';
+COLL_KCUT = 1.75;
+
+K_Ti = 6.96;
+DT = 1e-2;
+TMAX = 100;
+kymin = 0.05;
+Ny = 40;
+SIMID = 'linear_CBC_nu+PJ_scan_kT_6.96_SGGK'; % Name of the simulation
% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
RUN = 1;
-figure
-for P = [6 8 10]
+g_ky = zeros(numel(NU_a),numel(P_a),Ny/2);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+
+j = 1;
+for P = P_a
-i=1;
+i = 1;
for NU = NU_a
-
-%Set Up parameters
-for j = 1
+ %Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
TAU = 1.0; % e/i temperature ratio
- K_N = 2.22; K_Ne = K_N;
- K_Te = K_T; % Temperature '''
+ K_Ni = 2.22; K_Ne = K_Ni;
+ K_Te = K_Ti; % Temperature '''
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0e-1; % electron plasma beta
J = P/2;
+% J = 2;
PMAXE = P; JMAXE = J;
PMAXI = P; JMAXI = J;
- NX = 12; % real space x-gridpoints
- NY = 2; % '' y-gridpoints
+ NX = 2; % real space x-gridpoints
+ NY = Ny; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
- LY = 2*pi/ky_;
+ LY = 2*pi/kymin;
NZ = 16; % number of perpendicular planes (parallel grid)
- NPOL = 1; SG = 0;
+ NPOL = 2; SG = 0;
GEOMETRY= 's-alpha';
Q0 = 1.4; % safety factor
- SHEAR = 0.8; % magnetic shear (Not implemented yet)
- EPS = 0.18; % inverse aspect ratio
+ SHEAR = 0.8; % magnetic shear
+ KAPPA = 0.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+ EPS = 0.18; % inverse aspect ratio
SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
@@ -62,56 +76,84 @@ for j = 1
INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; N_HD = 4;
- MU_Z = 2.0; MU_P = 0.0; %
+ MU_Z = 0.2; MU_P = 0.0; %
MU_J = 0.0; LAMBDAD = 0.0;
- NOISE0 = 0.0e-5; % Init noise amplitude
- BCKGD0 = 1.0; % Init background
-GRADB = 1.0;CURVB = 1.0;
-end
-%%-------------------------------------------------------------------------
-% RUN
-setup
-if RUN
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
-end
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ GRADB = 1.0;CURVB = 1.0;
+ %%-------------------------------------------------------------------------
+ % RUN
+ setup
+ if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 2 3 1 0; cd ../../../wk'])
+ end
+
+ % Load results
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+
+ % linear growth rate (adapted for 2D zpinch and fluxtube)
+ options.TRANGE = [0.5 1]*data.Ts3D(end);
+ options.NPLOTS = 0; % 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);
-% Load results
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
-data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
-
-%linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 0;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
-[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);
-
- g_max(i) = gmax;
- k_max(i) = kmax;
+ g_ky(i,j,:) = lg.avg_g;
- g_avg(i) = lg.avg_g;
- g_std(i) = lg.std_g;
+ g_avg(i,j,:) = lg.avg_g;
+ g_std(i,j,:) = lg.std_g;
i = i + 1;
end
-%%
+j = j + 1;
+end
+
+if 1
+%% Study of the peak growth rate
+figure
-% plot(KT_a,max(g_max,0));
y_ = g_avg;
e_ = g_std;
+% filter to noisy data
y_ = y_.*(y_-e_>0);
e_ = e_ .* (y_>0);
-errorbar(NU_a,y_,e_,...
- 'LineWidth',1.2,...
- 'DisplayName',['(',num2str(P),',',num2str(P/2),')']);
-hold on;
-title(['$\kappa_T=$',num2str(K_T),' $k_y=$',num2str(ky_),' (CLOS = 0)']);
-legend('show'); xlabel('$\nu_{DGGK}$'); ylabel('$\gamma$');
-drawnow
+
+[y_,idx_] = max(g_avg,[],3);
+for i = 1:numel(idx_)
+ e_ = g_std(:,:,idx_(i));
+end
+
+subplot(121)
+for i = 1:numel(NU_a)
+ errorbar(P_a,y_(i,:),e_(i,:),...
+ 'LineWidth',1.2,...
+ 'DisplayName',['$\nu=$',num2str(NU_a(i))]);
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
+legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
+
+subplot(122)
+for j = 1:numel(P_a)
+ errorbar(NU_a,y_(:,j),e_(:,j),...
+ 'LineWidth',1.2,...
+ 'DisplayName',['(',num2str(P_a(j)),',',num2str(P_a(j)/2),')']);
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=$',num2str(data.ky(idx_))]);
+legend('show'); xlabel(['$\nu_{',CO,'}$']); ylabel('$\gamma$');
end
+if 0
+%% Pcolor of the peak
+figure;
+[XX_,YY_] = meshgrid(NU_a,P_a);
+pclr=pcolor(XX_,YY_,y_'); set(pclr,'EdgeColor','none');
+end
\ No newline at end of file
diff --git a/wk/Zpinch_coll_scan_kN_1.7.m b/wk/Zpinch_coll_scan_kN_1.7.m
index ad00ba14..4f470598 100644
--- a/wk/Zpinch_coll_scan_kN_1.7.m
+++ b/wk/Zpinch_coll_scan_kN_1.7.m
@@ -59,25 +59,144 @@ xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
set(gca,'Yscale','log');
end
+%% KN SCAN (4,2) 200x64, muHD = 0.1 to 1.0, NL = -1, DT=1e-2 nu = 0.1
if 0
%%
figure
nu = 0.1;
-% FCGK 4,2
+
+% SGGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40];
+Gavg_a = [0.0264 0.7829 3.4742 26.3825 73.3567];
+Gstd_a = [0.0039 0.1203 0.3892 4.9250 11.5890];
+errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Sugama (4,2)'); hold on
+
+% DGGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40];
+Gavg_a = [0.0080 0.0843 1.6752 28.0470 69.5804];
+Gstd_a = [0.0010 0.1303 0.3329 6.0057 11.4145];
+errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Dougherty (4,2)'); hold on
+
+% LRGK 4,2
kn_a = [1.60 1.80 2.00 2.20 2.40];
-Gavg_a = [1.11e-1 6.86e-1 3.44e-0 1.12e+1 2.87e+1];
-Gstd_a = [7.98e-3 1.10e-1 4.03e-1 2.03e+0 7.36e+0];
+Gavg_a = [0.8244 2.7805 6.4115 21.1632 62.2215];
+Gstd_a = [0.0621 0.2755 0.6950 3.6594 15.0836];
+errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Lorentz (4,2)'); hold on
+% FCGK 4,2
+kn_a = [1.60 1.80 2.00 2.20];
+Gavg_a = [0.5071 3.1178 10.1663 29.3928];
+Gstd_a = [0.0491 0.3237 1.1579 4.0400];
errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Coulomb (4,2)'); hold on
+% LDGKii 4,2
+% kn_a = [1.60 1.80 2.00 2.20 2.40];
+% Gavg_a = [1.11e-1 6.86e-1 3.44e-0 1.12e+1 2.87e+1];
+% Gstd_a = [7.98e-3 1.10e-1 4.03e-1 2.03e+0 7.36e+0];
+% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Landauii (4,2)'); hold on
+
% % Collisionless
% plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
%
+title('Transport w.r.t. gradient level at $\nu = 0.1$');
xlim([1.6 2.5]);
legend('show');
-xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
+xlabel('$L_n/R$'); ylabel('$\Gamma_x^\infty R/L_n$');
+end
+%% KN SCAN (4,2) 200x64, muHD = 0.1 to 1.0, NL = -1, DT=1e-2 nu = 0.01
+if 0
+%%
+figure
+Mksz = 10; Lnwd = 2.0;
+
+nu = 0.01;
+
+pnlty = 1;
+% SGGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
+Gavg_a = [0.0104 0.1735 0.6430 3.8566 37.0602 pnlty*0.0433 pnlty*0.5181 pnlty*1.7723 pnlty*14.5713 pnlty*85.4020];
+Gstd_a = [0.0045 0.0868 0.3587 1.8440 8.6218 pnlty*0.0214 pnlty*0.2620 pnlty*0.9196 pnlty* 4.1692 pnlty*21.5764];
+[~,Idx] = sort(kn_a);
+kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
+% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Sugama (4,2)'); hold on
+loglog(kn_a, Gavg_a./kn_a,'s','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Sugama (4,2)'); hold on
+
+% DGGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
+Gavg_a = [0.0010 0.1245 0.5219 1.7153 38.6814 pnlty*0.0020 pnlty*0.3208 pnlty*1.0034 pnlty*2.6310 pnlty*79.4978];
+Gstd_a = [0.0002 0.1154 0.5585 0.9062 12.814 pnlty*0.0005 pnlty*0.1853 pnlty*0.5899 pnlty*1.5473 pnlty*20.7551];
+[~,Idx] = sort(kn_a);
+kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
+% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Dougherty (4,2)'); hold on
+loglog(kn_a, Gavg_a./kn_a,'v','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Dougherty (4,2)'); hold on
+
+% LRGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
+Gavg_a = [0.1422 0.5389 1.2577 4.4465 50.4417 pnlty*0.2009 pnlty*0.6764 pnlty*1.8266 pnlty*9.0066 pnlty*83.4325];
+Gstd_a = [0.0119 0.1684 0.4304 2.8511 4.7634 pnlty*0.0434 pnlty*0.2677 pnlty*0.9563 pnlty*6.9812 pnlty*27.4134];
+[~,Idx] = sort(kn_a);
+kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
+% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Lorentz (4,2)'); hold on
+loglog(kn_a, Gavg_a./kn_a,'^','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Lorentz (4,2)'); hold on
+
+plot(1,1);
+% FCGK 4,2
+kn_a = [1.60 1.80 2.00 2.20 2.40 1.70 1.90 2.10 2.30 2.50];
+Gavg_a = [0.1914 0.8660 2.2324 5.3947 44.2023 pnlty*0.3803 pnlty*1.3516 pnlty*3.5982 pnlty*16.6212 pnlty*88.5276];
+Gstd_a = [0.0097 0.1009 0.4804 1.2491 9.8378 pnlty*0.0202 pnlty*0.3688 pnlty*0.8374 pnlty*3.1974 pnlty*19.0415];
+[~,Idx] = sort(kn_a);
+kn_a = kn_a(Idx); Gavg_a = Gavg_a(Idx); Gstd_a = Gstd_a(Idx);
+% errorbar(kn_a, Gavg_a./kn_a, Gstd_a./kn_a,'DisplayName','Coulomb (4,2)'); hold on
+loglog(kn_a, Gavg_a./kn_a,'d','MarkerSize',Mksz,'LineWidth',Lnwd,'DisplayName','Coulomb (4,2)'); hold on
+
+
+%
+% title('Transport w.r.t. gradient level at $\nu = 0.01$');
+xlim([1.55 2.55]);
+% legend('show');
+xlabel('$\kappa_N$'); ylabel('$\Gamma_x^\infty /\kappa_n$');
+grid on
+end
+
+if 0
+ %% Mixing length argument
+ figure;
+ nu = 0.0;
+ % collisionless
+ kn= [1.6 2.0 2.5 ];
+ k = [0.65 0.60 0.50];
+ g = [0.19 0.44 0.85];
+ y = g.^2./k.^3;
+ plot(kn,y*(9.6/y(3)),'--k'); hold on;
+
+ nu = 0.01;
+ % DG
+ kn= [1.6 2.0 2.5 ];
+ k = [0.38 0.60 0.66];
+ g = [0.13 0.50 0.90];
+ y = g.^2./k.^3;
+ plot(kn,y*(9.6/y(3)),'--r');
+ % SG
+ kn= [1.6 2.0 2.5 ];
+ k = [0.41 0.63 0.69];
+ g = [0.16 0.51 0.89];
+ y = g.^2./k.^3;
+ plot(kn,y*(9.6/y(3)),'--b');
+ % LR
+ kn= [1.6 2.0 2.5 ];
+ k = [0.41 0.57 0.60];
+ g = [0.18 0.50 0.87];
+ y = g.^2./k.^3;
+ plot(kn,y*(9.6/y(3)),'--y');
+ % LD
+ kn= [1.6 2.0 2.5 ];
+ k = [0.38 0.53 0.57];
+ g = [0.15 0.47 0.85];
+ y = g.^2./k.^3;
+ plot(kn,y*(9.6/y(3)),'--g');
+
end
\ No newline at end of file
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 53f11c4f..a35d30f9 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -1,8 +1,8 @@
-helazdir = '/home/ahoffman/HeLaZ/';
-addpath(genpath([helazdir,'matlab'])) % ... add
-addpath(genpath([helazdir,'matlab/plot'])) % ... add
-addpath(genpath([helazdir,'matlab/compute'])) % ... add
-addpath(genpath([helazdir,'matlab/load'])) % ... add
+gyacomodir = '/home/ahoffman/gyacomo/';
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
@@ -12,14 +12,17 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.5/';
% 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_2b_mu_5e-2/';
+% folder = '/misc/gene_results/Z-pinch/HP_fig_2a_mu_1e-2/';
+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/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_00/';
+% folder = '/misc/gene_results/Z-pinch/kN_2.5_HD_transport_spectrum_00/';
+
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/196x96x20x32x16_02/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
@@ -39,6 +42,7 @@ options.TAVG_1 = gene_data.Ts3D(end); % Averaging times duration
options.NMVA = 1; % Moving average for time traces
options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x, Q_x)
options.INTERP = 1;
+options.NCUT = 4; % Number of cuts for averaging and error estimation
gene_data.FIGDIR = folder;
fig = plot_radial_transport_and_spacetime(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -57,15 +61,15 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
% options.NAME = 'Q_x';
-options.NAME = '\phi';
-% options.NAME = 'n_i';
+% options.NAME = '\phi';
+options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [325 400];
+options.TIME = [1:10];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -76,12 +80,11 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'v_y';
-options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -118,7 +121,7 @@ end
if 0
%% Show f_i(vpar,mu)
-options.times = 250:500;
+options.times = 600:5000;
options.specie = 'i';
options.PLT_FCT = 'contour';
options.folder = folder;
@@ -131,7 +134,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [4000 8000];
+options.TIME = [60 10000];
options.NORM =1;
% options.NAME = '\phi';
% options.NAME = 'n_i';
@@ -150,7 +153,7 @@ if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 1;
-options.TIME = 100:700;
+options.TIME = 1:700;
options.NMA = 1;
options.NMODES = 15;
options.iz = 'avg';
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 2a875f63..8d494f0a 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -1,5 +1,5 @@
%% UNCOMMENT FOR TUTORIAL
-% gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
+gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
% resdir = '.'; %Name of the directory where the results are located
% JOBNUMMIN = 00; JOBNUMMAX = 10;
%%
@@ -31,7 +31,7 @@ if 1
i_ = 1;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
-options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
+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
@@ -40,46 +40,49 @@ options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag
options.INTERP = 0;
options.RESOLUTION = 256;
fig = plot_radial_transport_and_spacetime(data,options);
-save_figure(data,fig,'.png')
+% save_figure(data,fig,'.png')
end
if 0
%% statistical transport averaging
-for i_ = 1:2:21
+gavg =[]; gstd = [];
+for i_ = 3: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))])
options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)];
options.NPLOTS = 0;
-fig = statistical_transport_averaging(data,options);
+[fig, gavg_, gstd_] = statistical_transport_averaging(data,options);
+gavg = [gavg gavg_]; gstd = [gstd gstd_];
end
+disp(gavg); disp(gstd);
end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = '\omega_z';
% options.NAME = 'N_i^{00}';
-% 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.NAME = 'n_i';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [000:0.1:7000];
+options.TIME = [000:1:8000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
create_film(data,options,'.gif')
end
-if 1
+if 0
%% 2D snapshots
% Options
options.INTERP = 1;
@@ -96,7 +99,7 @@ options.PLAN = 'xy';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [1000 1800 2500 3000 4000];
+options.TIME = [1000 3000 5000 8000 10000];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
@@ -119,16 +122,16 @@ end
if 0
%% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
-options.SPAR = linspace(-3,3,32)+(6/127/2);
-options.XPERP = linspace( 0,6,32);
-% options.SPAR = gene_data.vp';
-% options.XPERP = gene_data.mu';
+% options.SPAR = linspace(-3,3,32)+(6/127/2);
+% options.XPERP = linspace( 0,6,32);
+options.SPAR = gene_data.vp';
+options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [250 600];
-options.PLT_FCT = 'pcolor';
-options.ONED = 0;
+options.T = [1500 5000];
+options.PLT_FCT = 'contour';
+options.ONED = 1;
options.non_adiab = 0;
-options.SPECIE = 'i';
+options.SPECIE = 'e';
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')
@@ -175,8 +178,8 @@ options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME ='\Gamma_x';
% options.NAME ='s_y';
-options.COMPX = 'avg';
-options.COMPY = 'avg';
+options.COMPX = 1;
+options.COMPY = 2;
options.COMPZ = 1;
options.COMPT = 1;
options.MOVMT = 1;
@@ -220,11 +223,11 @@ end
if 0
%% linear growth rate for 3D Zpinch
-trange = [5 15];
+trange = [5 20];
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,'.png')
+% save_figure(data,fig,'.png')
end
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index cdd4074d..11cfa44a 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';
@@ -197,16 +197,18 @@ resdir ='';
% % resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
% resdir = 'Zpinch_rerun/kN_2.2_SGGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
-%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0)
+
+%% Convergence cases kN = (1.6 2.2) nu = (0.01 1.0), SGGK
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x7x4_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x9x5_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x11x6_nu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x3x2_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x5x3_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_mu_0.01';
-resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
+% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x3x2_nu_0.1';
% resdir = 'Zpinch_rerun/convcoll_Kn_1.6_200x32x5x3_nu_0.1';
@@ -218,6 +220,22 @@ resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_mu_0.01';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x7x4_nu_0.1';
% resdir = 'Zpinch_rerun/convcoll_Kn_2.2_200x32x9x5_nu_0.1';
+%% kN scans with nu=0.1 and nu=0.01
+% resdir = 'Zpinch_rerun/SGGK_kN_scan_200x64x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/DGGK_kN_scan_200x64x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/LRGK_kN_scan_200x64x5x3_nu_0.1';
+% resdir = 'Zpinch_rerun/LDGK_kN_scan_200x64x5x3_nu_0.1';
+%
+% resdir = 'Zpinch_rerun/SGGK_kN_scan_200x64x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/DGGK_kN_scan_200x64x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/even_DGGK_kN_scan_200x64x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/LRGK_kN_scan_200x64x5x3_nu_0.01';
+% resdir = 'Zpinch_rerun/LDGK_kN_scan_200x64x5x3_nu_0.01';
+
+% 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/LDGK_kN_1.9_200x64x5x3_nu_0.01';
JOBNUMMIN = 00; JOBNUMMAX = 03;
resdir = ['results/',resdir];
diff --git a/wk/lin_3D_Zpinch.m b/wk/lin_3D_Zpinch.m
new file mode 100644
index 00000000..a1de410d
--- /dev/null
+++ b/wk/lin_3D_Zpinch.m
@@ -0,0 +1,194 @@
+%% QUICK RUN SCRIPT
+% This script create a directory in /results and run a simulation directly
+% 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_3D_Zpinch'; % Name of the simulation
+SIMID = 'NL_3D_zpinch'; % 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
+PROGDIR = '/home/ahoffman/gyacomo/';
+% EXECNAME = 'gyacomo_1.0';
+EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.05; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 2.0; % ele Density '''
+K_Te = K_Ne/4; % ele Temperature '''
+K_Ni = K_Ne; % ion Density gradient drive
+K_Ti = K_Ni/4; % ion Temperature '''
+% SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.005; % electron plasma beta
+%% GRID PARAMETERS
+P = 2;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 128; % real space x-gridpoints
+NY = 32; % '' y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain
+LY = 2*pi/0.1; % Size of the squared frequency domain
+NZ = 20; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+GEOMETRY= 'Z-pinch';
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+KAPPA = 0.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 20; % 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 = 5; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'DG';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax 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
+%% OUTPUTS
+W_DOUBLE = 1;
+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.0; % Hyper diffusivity cutoff ratio
+MU = 0.0; % Hyperdiffusivity coefficient
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+MU_X = MU; %
+MU_Y = MU; %
+N_HD = 4;
+MU_Z = 2.0; %
+MU_P = 0.0; %
+MU_J = 0.0; %
+LAMBDAD = 0.0;
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
+GRADB = 1.0;
+CURVB = 1.0;
+COLL_KCUT = 1000;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+if NZ > 1
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 0; cd ../../../wk'])
+else
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+end
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 1
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+options.TRANGE = [0.5 1]*data.Ts3D(end);
+options.NPLOTS = 3; % 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 1
+%% Ballooning plot
+options.time_2_plot = [120];
+options.kymodes = [0.3];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 1;
+options.ST = 1;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
+if 0
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+options.NMA = 1;
+options.NMODES = 1;
+options.iz = 'avg';
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 1
+%% RH TEST
+ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
+end
diff --git a/wk/lin_EPY.m b/wk/lin_ETPY.m
similarity index 87%
rename from wk/lin_EPY.m
rename to wk/lin_ETPY.m
index 2ddd756a..101f81c3 100644
--- a/wk/lin_EPY.m
+++ b/wk/lin_ETPY.m
@@ -1,9 +1,9 @@
-%% QUICK RUN SCRIPT for linear entropy mode in a Zpinch
+%% QUICK RUN SCRIPT for linear entropy mode (ETPY) in a Zpinch
% This script create a directory in /results and run a simulation directly
% 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_EPY'; % Name of the simulation
+SIMID = 'lin_ETPY'; % Name of the simulation
% SIMID = 'dbg'; % Name of the simulation
RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
@@ -16,9 +16,9 @@ EXECNAME = 'gyacomo';
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_Ne = 1.6; % ele Density '''
+K_Ne = 2.5; % ele Density '''
K_Te = K_Ne/4; % ele Temperature '''
K_Ni = K_Ne; % ion Density gradient drive
K_Ti = K_Ni/4; % ion Temperature '''
@@ -26,7 +26,7 @@ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 2;
+P = 4;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
@@ -35,7 +35,7 @@ JMAXI = J; % "
NX = 2; % real space x-gridpoints
NY = 100; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 120;%2*pi/0.05; % Size of the squared frequency domain
+LY = 200;%2*pi/0.05; % Size of the squared frequency domain
NZ = 1; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
@@ -43,14 +43,17 @@ SG = 0; % Staggered z grids option
GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
+KAPPA = 0.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
TMAX = 500; % Maximal time unit
-DT = 1e-2; % Time step
+DT = 5e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
-SPS3D = 1/5; % Sampling per time unit for 2D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
SPS5D = 1/5; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
@@ -58,7 +61,7 @@ JOB2LOAD= -1;
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
% Collision operator
% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-CO = 'SG';
+CO = 'LD';
GKCO = 1; % gyrokinetic operator
ABCO = 1; % interspecies collisions
INIT_ZF = 0; ZF_AMP = 0.0;
@@ -85,18 +88,19 @@ MU_Z = 2.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
-NOISE0 = 0.0e-5; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
+NOISE0 = 1.0e-5; % Init noise amplitude
+BCKGD0 = 0.0; % Init background
GRADB = 1.0;
CURVB = 1.0;
+COLL_KCUT = 1000;
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',PROGDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
@@ -157,7 +161,7 @@ options.kzky = 0;
[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
save_figure(data,fig)
end
-if 1
+if 0
%% Mode evolution
options.NORMALIZED = 0;
options.K2PLOT = 10;
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 3ebe622d..6f18bee5 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -7,7 +7,7 @@ SIMID = 'lin_ITG'; % Name of the simulation
RUN = 0; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
-HELAZDIR = '/home/ahoffman/gyacomo/';
+gyacomodir = '/home/ahoffman/gyacomo/';
% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
@@ -103,14 +103,14 @@ 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,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 1 2 3 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'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+LOCALDIR = [gyacomodir,'results/',filename,'/'];
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 00;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -118,9 +118,10 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
%% Short analysis
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
-trange = [0.5 1]*data.Ts3D(end);
-nplots = 3;
-lg = compute_fluxtube_growth_rate(data,trange,nplots);
+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);
--
GitLab