diff --git a/matlab/compute/invert_kxky_to_kykx_gene_results.m b/matlab/compute/invert_kxky_to_kykx_gene_results.m
index 904bc6d028d346bd099d2dcafa29f57e74599dfc..462846b7a85aa00d4f964901258455899898df68 100644
--- a/matlab/compute/invert_kxky_to_kykx_gene_results.m
+++ b/matlab/compute/invert_kxky_to_kykx_gene_results.m
@@ -12,7 +12,13 @@ data.TEMP_I = rotate(data.TEMP_I);
data.Ny = data.Nky*2-1;
data.Nx = data.Nkx;
-dkx = data.kx(2); dky = data.ky(2);
+if numel(data.kx)>1
+ dkx = data.kx(2);
+else
+ dkx = 1;
+end
+
+dky = data.ky(2);
Lx = 2*pi/dkx; Ly = 2*pi/dky;
x = linspace(-Lx/2,Lx/2,data.Nx+1); x = x(1:end-1);
y = linspace(-Ly/2,Ly/2,data.Ny+1); y = y(1:end-1);
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 108efd9617cc453f3c9a23d1d187972328930838..dda3387ecda1f1598d7d0b95dc52f4fc88ccc3aa 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -17,6 +17,9 @@ end
FRAMES = unique(FRAMES);
t = DATA.Ts3D(FRAMES);
+% time window where we measure the growth
+it1 = floor(numel(t)/2);
+it2 = numel(t);
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
@@ -77,7 +80,8 @@ for i = 1:2
amp = MODES;
i_=1;
for ik = MODES
- gr = polyfit(t,log(plt(field,ik)),1);
+ to_measure = log(plt(field,ik));
+ gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
gamma(i_) = gr(1);
amp(i_) = gr(2);
i_=i_+1;
@@ -101,14 +105,18 @@ for i = 1:2
if MODES_SELECTOR == 2
semilogy(t,plt(field,ikzf),'--k');
end
+ %plot the time window where the gr are measured
+ plot(t(it1)*[1 1],[1e-10 1],'--k')
+ plot(t(it2)*[1 1],[1e-10 1],'--k')
xlim([t(1) t(end)]); %ylim([1e-5 1])
xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
title('Measure time window')
subplot(2,3,3+3*(i-1))
- plot(k(MODES),gamma); hold on;
+ plot(k(MODES),gamma,...
+ 'DisplayName',['(',num2str(DATA.Pmaxi-1),',',num2str(DATA.Jmaxi-1),')']); hold on;
for i_ = 1:numel(mod2plot)
- plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:))
+ plot(k(MODES(mod2plot(i_))),gamma(mod2plot(i_)),'x','color',clr_(i_,:));
end
if MODES_SELECTOR == 2
plot(k(ikzf),gamma(ikzf),'ok');
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 1ccdb54bc63ef4761d18168ba12030238fb608c3..641e6dd9aff55ea85d377d20fc3c1e984c3a9e22 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -21,7 +21,12 @@ DATA.Ny = DATA.Nky*2-1;
DATA.z = h5read([folder,coofile],'/coord/z');
DATA.Nz = numel(DATA.z);
-dkx = DATA.kx(2); dky = DATA.ky(2);
+if numel(DATA.kx)>1
+ dkx = DATA.kx(2);
+else
+ dkx = 1;
+end
+dky = DATA.ky(2);
Lx = 2*pi/dkx; Ly = 2*pi/dky;
x = linspace(-Lx/2,Lx/2,DATA.Nx+1); x = x(1:end-1);
y = linspace(-Ly/2,Ly/2,DATA.Ny+1); y = y(1:end-1);
@@ -48,14 +53,14 @@ for jt = 1:numel(DATA.Ts3D)
t = DATA.Ts3D(jt);
[~, it] = min(abs(DATA.Ts3D-t));
%
- tmp = h5read([folder,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
- DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-%
- tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
- DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
-%
- 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,momfile],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+% DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+% %
+% tmp = h5read([folder,momfile],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
+% DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+% %
+% 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/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 8fb67265b2114225e6482cc0fa57574bf97c1c93..55699d4e9a9edeef1da545e39173fedb0e40a007 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -41,7 +41,7 @@
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
path = '/home/ahoffman/gene/linear_CBC_results/';
-fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
+% fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
% fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt';
@@ -52,8 +52,38 @@ fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
% fname = 'CBC_KT_11_20x1x16x24x10_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_11_20x1x32x30x14_Lv_3_Lw_12.txt';
% fname = 'CBC_ky_0.3_20x1x16x24x10_Lv_3_Lw_12_nuv_1e-3.txt';
+%----------Shearless CBC
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_kine.txt';
+% fname = 'CBC_salpha_s0_nz_24_nv_48_nw_16_kine_beta_1e-4.txt';
+% fname = 'CBC_miller_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'CBC_miller_s0_nz_24_nv_48_nw_16_kine.txt';
+%----------Shearless pITG
+% fname = 'pITG_salpha_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_48_nw_16_adiabe.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_48_nw_16_kine.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_48_nw_16_kine.txt';
+%----------Convergence nvpar shearless pITG
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_16_kine.txt';
+% fname = 'pITG_salpha_s0_nz_24_nv_scan_nw_24_adiabe.txt';
+% fname = 'pITG_miller_s0_nz_24_nv_scan_nw_24_adiabe.txt';
+%----------Convergence nvpar shearless CBC
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'CBC_miller_nz_24_nv_scan_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_nz_24_nv_scan_nw_16_kine.txt';
+fname = 'CBC_miller_nz_24_nv_scan_nw_16_kine.txt';
+%---------- CBC
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+% fname = 'CBC_salpha_nx_8_nz_24_nv_36_nw_16_kine.txt';
+% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_adiabe.txt';
+% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_adiabe.txt';
+% fname = 'CBC_miller_nx_20_nz_32_nv_32_nw_12_kine.txt';
+% fname = 'CBC_miller_nx_8_nz_24_nv_36_nw_16_kine.txt';
data_ = load([path,fname]);
figure
-plot(data_(:,2),data_(:,3),'-dk','DisplayName','GENE'); hold on;
-plot(data_(:,2),data_(:,4),'--*k','DisplayName','GENE');
\ No newline at end of file
+plot(data_(:,2),data_(:,3),'-dk','DisplayName',fname); hold on;
+plot(data_(:,2),data_(:,4),'--*k','DisplayName',fname);
\ No newline at end of file
diff --git a/matlab/plot/plot_fa.m b/matlab/plot/plot_fa.m
index efdbd15f0f86bde554499b9b996fb1506c9b958a..2395489c00ae1b8ddc38d69c99c60537ff193ae9 100644
--- a/matlab/plot/plot_fa.m
+++ b/matlab/plot/plot_fa.m
@@ -26,17 +26,27 @@ if OPTIONS.ONED
end
else
- [SS,XX,FFa] = compute_fa_2D(DATA, OPTIONS); sz = size(SS);
+ [SS,XX,FFa] = compute_fa_2D(DATA, OPTIONS); sz = size(SS); FFa = FFa';
[~,it] = min(abs(OPTIONS.T-DATA.Ts5D));
switch OPTIONS.PLT_FCT
case 'contour'
- contour(SS,XX,FFa',sum(sz)/2);
+ contour(SS,XX,FFa,sum(sz)/2);
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
case 'pcolor'
- pclr = pcolor(SS,XX,FFa'); set(pclr, 'edgecolor','none'); shading interp
+ pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
case 'contourf'
- contourf(SS,XX,FFa',sum(sz)/2);
+ contourf(SS,XX,FFa,sum(sz)/2);
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
+ case 'surf'
+ surf(SS,XX,FFa);
+ xlabel('$v_\parallel$'); ylabel('$\mu$');
+ case 'surfvv'
+ surf([SS(end:-1:1,:) SS ],[-sqrt(XX(end:-1:1,:)) sqrt(XX)],[FFa(end:-1:1,:) FFa]);
+ xlabel('$v_\parallel$'); ylabel('$v_\perp$');
+ xlim([min(OPTIONS.SPAR) max(OPTIONS.SPAR)]);
+ ylim(sqrt(max(OPTIONS.XPERP))*[-1 1]);
end
- xlabel('$v_\parallel$'); ylabel('$\mu$');
legend(['$\langle |f_',OPTIONS.SPECIE,'|^2\rangle_{xy}^{1/2}$',zcomp])
if numel(it) == 1
title(['HeLaZ''$\langle |f_',OPTIONS.SPECIE...
diff --git a/matlab/plot/plot_fa_gene.m b/matlab/plot/plot_fa_gene.m
index 9792d557cf8ce6de76592722717b33a63da4c614..9743d80a3615d5b8e2a38239985c97851069b588 100644
--- a/matlab/plot/plot_fa_gene.m
+++ b/matlab/plot/plot_fa_gene.m
@@ -73,6 +73,10 @@ switch specie
pclr = contourf(SS,XX,FFa);
case 'pcolor'
pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
+ case 'surf'
+ surf(SS,XX,FFa);
+ case 'surfvv'
+ surf([SS; SS(:,end-1:-1:1)],sqrt([XX; XX(:,end-1:-1:1)]),[FFa; FFa(:,end-1:-1:1)]);
end
case 'i'
name = '$f_i(v_\parallel,\mu_p)$';
@@ -85,6 +89,13 @@ switch specie
pclr = contourf(SS,XX,FFa,(nvp+nmu)/2);
case 'pcolor'
pclr = pcolor(SS,XX,FFa); set(pclr, 'edgecolor','none'); shading interp
+ case 'surf'
+ surf(SS,XX,FFa);
+ case 'surfvv'
+ surf([SS(:,end:-1:1) SS ],[-sqrt(XX(:,end:-1:1)) sqrt(XX)],[FFa(:,end:-1:1) FFa]);
+ xlabel('$v_\parallel$'); ylabel('$v_\perp$');
+ xlim([min(vp) max(vp)]);
+ ylim(sqrt(max(mu))*[-1 1]);
end
end
if numel(TIMES) == 1
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 377c80fcaaf213f3867d70b5d423eaa133fedb9d..0eb8c09645229df907cf511aaa28742d19fa48e4 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,7 +1,10 @@
function [ fig ] = plot_metric( data, options )
-names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
- 'gxx','gxy','gyy','gyz','gzz','hatB','hatR','hatZ'};
+% names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
+% 'gxx','gxy','gxz','gyy','gyz','gzz','hatB','hatR','hatZ'};
+names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
+ 'hatB', 'gradxB', 'gradyB', 'gradzB',...
+ 'Jacobian','hatR','hatZ','gradz_coeff'};
geo_arrays = zeros(2,data.Nz,numel(names));
for i_ = 1:numel(names)
@@ -13,19 +16,19 @@ if NPLOT > 0
fig = figure;
if options.SHOW_METRICS
subplot(311)
- for i = 1:5
+ for i = 1:6
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
end
xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
subplot(312)
- for i = 6:10
+ for i = 7:10
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
end
xlim([min(data.z),max(data.z)]); legend('show');
subplot(313)
- for i = 11:13
+ for i = 11:14
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
end
xlim([min(data.z),max(data.z)]); legend('show');
diff --git a/matlab/setup.m b/matlab/setup.m
index 8e6245b7e6ed9205d69c7630ff6c3b0d01dc8442..8f3981190fb4f2e33e45110824b31dcafc951c93 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -136,7 +136,11 @@ if (NZ > 1) %3D case
end
switch GEOMETRY
case 'circular'
- geo_ = [geo_,'_circ_'];
+ geo_ = [geo_,'_cir_'];
+ case 's-alpha'
+ geo_ = [geo_,'_s-a_'];
+ case 'miller'
+ geo_ = [geo_,'_mil_'];
end
% put everything together in the param character chain
@@ -175,9 +179,9 @@ end
if ~exist(BASIC.RESDIR, 'dir')
mkdir(BASIC.RESDIR)
end
-if ~exist(BASIC.MISCDIR, 'dir')
-mkdir(BASIC.MISCDIR)
-end
+% if ~exist(BASIC.MISCDIR, 'dir')
+% mkdir(BASIC.MISCDIR)
+% end
%% Compile and WRITE input file
INPUT = write_fort90(OUTPUTS,GRID,GEOM,MODEL,COLL,INITIAL,TIME_INTEGRATION,BASIC);
nproc = 1;
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
index 9a07ac6e5d6b8359fda849fcc37634b31c10ebae..84b33df4ab7446e3479d46d057af4240e64dedf0 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_CO_scan.m
@@ -5,94 +5,129 @@ addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_dbg';
EXECNAME = 'gyacomo';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
-NU_a = [0.005 0.01:0.01:0.1];
-% P_a = [2 4 6 8 10 12 16];
-% NU_a = 0.1;
-P_a = [2 4 8 10 12 16 20];
-
-
-CO = 'SG';
+% SIMID = 'linear_CBC_nu+PJ_scan_kT_6.96_SGGK'; % Name of the simulation
+SIMID = 'convergence_pITG_dbg'; % Name of the simulation
+% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
+RERUN = 1; % If you want to rerun the sim (bypass the check of existing data)
+
+% NU_a = [0.0 0.01 0.02 0.05 0.1];
+NU_a = [0.0];
+P_a = [30];
+% P_a = [2 4:4:36];
+J_a = floor(P_a/2);
+% collision setting
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
COLL_KCUT = 1.75;
-
-K_Ti = 6.96;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_Ne = 0*2.22; % ele Density '''
+K_Te = 0*6.96; % ele Temperature '''
+K_Ni = 0*2.22; % ion Density gradient drive
+K_Ti = 6.96; % ion Temperature '''
+% Geometry
+GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
+SHEAR = 0.0; % magnetic shear
+% time and numerical grid
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 = 0;
-
-g_ky = zeros(numel(NU_a),numel(P_a),Ny/2);
+TMAX = 50;
+kymin = 0.4;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
g_avg= g_ky*0;
g_std= g_ky*0;
j = 1;
for P = P_a
-
i = 1;
for NU = NU_a
- %Set Up parameters
- CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
- TAU = 1.0; % e/i temperature ratio
- K_Ni = 2.22; K_Ne = K_Ni;
- K_Te = K_Ti; % Temperature '''
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % 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 = 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;
+ %% GRID PARAMETERS
+% P = 20;
+ J = floor(P/2);
+ PMAXE = P; % Hermite basis size of electrons
+ JMAXE = J; % Laguerre "
+ PMAXI = P; % " ions
+ JMAXI = J; % "
NX = 2; % real space x-gridpoints
- NY = Ny; % '' y-gridpoints
- LX = 2*pi/0.1; % Size of the squared frequency domain
- LY = 2*pi/kymin;
- NZ = 16; % number of perpendicular planes (parallel grid)
- NPOL = 2; SG = 0;
- GEOMETRY= 's-alpha';
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % 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';
+ EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
- SHEAR = 0.8; % magnetic shear
- KAPPA = 0.0; % elongation
+ KAPPA = 1.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 = -1; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
+ % PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+ PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ SPS0D = 1; % Sampling per time unit for 2D arrays
+ SPS2D = -1; % Sampling per time unit for 2D arrays
+ SPS3D = 1; % Sampling per time unit for 2D arrays
+ SPS5D = 1/2; % 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)
- GKCO = 1; % gyrokinetic operator
+ % Collision 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
+ INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% 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 = 0.2; MU_P = 0.0; %
- MU_J = 0.0; LAMBDAD = 0.0;
+ MU_Y = MU; %
+ N_HD = 4;
+ MU_Z = 0.2; %
+ 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;
- %%-------------------------------------------------------------------------
- % RUN
+ 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'])
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if (RERUN || isempty(data.NU_EVOL))
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
- % Load results
+ % Load results after trying to run
filename = [SIMID,'/',PARAMS,'/'];
LOCALDIR = [gyacomodir,'results/',filename,'/'];
data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
@@ -101,17 +136,29 @@ for NU = NU_a
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);
-
+% lg = compute_fluxtube_growth_rate(data,options);
+% [gmax, kmax] = max(lg.g_ky(:,end));
+% [gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+% g_ky(i,j,:) = g_ky;
+%
+% g_avg(i,j,:) = lg.avg_g;
+% g_std(i,j,:) = lg.std_g;
- g_ky(i,j,:) = lg.avg_g;
+ [~,it1] = min(abs(data.Ts3D-0.5*data.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data.Ts3D-1.0*data.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 1:NY/2+1
+ field = squeeze(sum(abs(data.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t);
+ gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
+ g_ky(i,j,ik) = gr(1);
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
- g_avg(i,j,:) = lg.avg_g;
- g_std(i,j,:) = lg.std_g;
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
+
i = i + 1;
end
@@ -122,23 +169,27 @@ if 1
%% Study of the peak growth rate
figure
-y_ = g_avg;
-e_ = g_std;
+y_ = g_ky;
+e_ = 0.05;
% filter to noisy data
y_ = y_.*(y_-e_>0);
e_ = e_ .* (y_>0);
-[y_,idx_] = max(g_avg,[],3);
+[y_,idx_] = max(g_ky,[],3);
for i = 1:numel(idx_)
e_ = g_std(:,:,idx_(i));
end
-colors_ = summer(numel(NU_a));
+colors_ = lines(numel(NU_a));
subplot(121)
for i = 1:numel(NU_a)
- errorbar(P_a,y_(i,:),e_(i,:),...
- 'LineWidth',1.2,...
+% errorbar(P_a,y_(i,:),e_(i,:),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
+% 'color',colors_(i,:));
+ plot(P_a,y_(i,:),'s-',...
+ 'LineWidth',2.0,...
'DisplayName',['$\nu=$',num2str(NU_a(i))],...
'color',colors_(i,:));
hold on;
@@ -149,9 +200,13 @@ legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
colors_ = jet(numel(P_a));
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),')'],...
+% errorbar(NU_a,y_(:,j),e_(:,j),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+% 'color',colors_(j,:));
+ plot(NU_a,y_(:,j),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
'color',colors_(j,:));
hold on;
end
diff --git a/wk/CBC_nu_CO_scan_miller.m b/wk/CBC_nu_CO_scan_miller.m
new file mode 100644
index 0000000000000000000000000000000000000000..87b104ed61b6fa3f67b370ba44845991492361a1
--- /dev/null
+++ b/wk/CBC_nu_CO_scan_miller.m
@@ -0,0 +1,227 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_dbg';
+% EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo_alpha';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+% SIMID = 'linear_CBC_nu+PJ_scan_kT_6.96_SGGK'; % Name of the simulation
+SIMID = 'convergence_pITG_miller'; % Name of the simulation
+% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
+RUN = 0; % to make sure it does not try to run
+RERUN = 0; % rerun if the data does not exist
+
+NU_a = [0.02];
+P_a = [24];
+% NU_a = [0.0 0.01 0.02 0.05 0.1];
+% P_a = [2 4:4:28];
+% P_a = [24:4:28];
+J_a = floor(P_a/2);
+% collision setting
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_Ne = 1*2.22; % ele Density '''
+K_Te = 1*6.96; % ele Temperature '''
+K_Ni = 1*2.22; % ion Density gradient drive
+K_Ti = 6.96; % ion Temperature '''
+% Geometry
+GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+% DT = 2e-3;
+DT = 5e-4;
+TMAX = 50;
+kymin = 0.4;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+
+j = 1;
+for P = P_a
+i = 1;
+for NU = NU_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % 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)
+ %% GRID PARAMETERS
+% P = 20;
+ J = floor(P/2);
+ PMAXE = P; % Hermite basis size of electrons
+ JMAXE = J; % Laguerre "
+ PMAXI = P; % " ions
+ JMAXI = J; % "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % 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';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ % PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+ PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ SPS0D = 1; % Sampling per time unit for 2D arrays
+ SPS2D = -1; % Sampling per time unit for 2D arrays
+ SPS3D = 1; % Sampling per time unit for 2D arrays
+ SPS5D = 1/2; % 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
+ 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= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% 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 = 0.2; %
+ 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;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if ((RERUN || isempty(data.NU_EVOL)) && RUN)
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
+ end
+
+ % Load results after trying to run
+ 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);
+% g_ky(i,j,:) = g_ky;
+%
+% g_avg(i,j,:) = lg.avg_g;
+% g_std(i,j,:) = lg.std_g;
+
+ [~,it1] = min(abs(data.Ts3D-0.5*data.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data.Ts3D-1.0*data.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 1:NY/2+1
+ field = squeeze(sum(abs(data.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t);
+ gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
+ g_ky(i,j,ik) = gr(1);
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
+
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 1
+%% Study of the peak growth rate
+figure
+
+y_ = g_ky;
+e_ = 0.05;
+
+% filter to noisy data
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+
+[y_,idx_] = max(g_ky,[],3);
+for i = 1:numel(idx_)
+ e_ = g_std(:,:,idx_(i));
+end
+
+colors_ = lines(numel(NU_a));
+subplot(121)
+for i = 1:numel(NU_a)
+% errorbar(P_a,y_(i,:),e_(i,:),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
+% 'color',colors_(i,:));
+ plot(P_a,y_(i,:),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
+ 'color',colors_(i,:));
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
+legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
+
+colors_ = jet(numel(P_a));
+subplot(122)
+for j = 1:numel(P_a)
+% errorbar(NU_a,y_(:,j),e_(:,j),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+% 'color',colors_(j,:));
+ plot(NU_a,y_(:,j),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+ 'color',colors_(j,:));
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
+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/CBC_nu_CO_scan_salpha.m b/wk/CBC_nu_CO_scan_salpha.m
new file mode 100644
index 0000000000000000000000000000000000000000..6283034931f640c46e348a05f881c51e00daa340
--- /dev/null
+++ b/wk/CBC_nu_CO_scan_salpha.m
@@ -0,0 +1,226 @@
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
+% EXECNAME = 'gyacomo_debug';
+% EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo_alpha';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'convergence_CBC_salpha'; % Name of the simulation
+% SIMID = 'dbg'; % Name of the simulation
+RUN = 0; % to make sure it does not try to run
+RERUN = 0; % rerun if the data does not exist
+
+% NU_a = [0.05];
+% P_a = [30];
+NU_a = [0.0 0.01 0.02 0.05 0.1];
+% P_a = [2 4:4:12];
+P_a = [2 4:4:28];
+J_a = floor(P_a/2);
+% collision setting
+CO = 'DG';
+GKCO = 0; % gyrokinetic operator
+COLL_KCUT = 1.75;
+% model
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_Ne = 1*2.22; % ele Density '''
+K_Te = 1*6.96; % ele Temperature '''
+K_Ni = 1*2.22; % ion Density gradient drive
+K_Ti = 6.96; % ion Temperature '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+% DT = 2e-3;
+DT = 5e-4;
+TMAX = 50;
+kymin = 0.4;
+NY = 2;
+% arrays for the result
+g_ky = zeros(numel(NU_a),numel(P_a),NY/2+1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+
+j = 1;
+for P = P_a
+i = 1;
+for NU = NU_a
+ %% PHYSICAL PARAMETERS
+ TAU = 1.0; % e/i temperature ratio
+ % 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)
+ %% GRID PARAMETERS
+% P = 20;
+ J = floor(P/2);
+ PMAXE = P; % Hermite basis size of electrons
+ JMAXE = J; % Laguerre "
+ PMAXI = P; % " ions
+ JMAXI = J; % "
+ NX = 8; % real space x-gridpoints
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/kymin; % Size of the squared frequency domain
+ NZ = 24; % 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';
+ EPS = 0.18; % inverse aspect ratio
+ Q0 = 1.4; % safety factor
+ KAPPA = 1.0; % elongation
+ DELTA = 0.0; % triangularity
+ ZETA = 0.0; % squareness
+ PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
+ SHIFT_Y = 0.0;
+ %% TIME PARMETERS
+ SPS0D = 1; % Sampling per time unit for 2D arrays
+ SPS2D = -1; % Sampling per time unit for 2D arrays
+ SPS3D = 1; % Sampling per time unit for 2D arrays
+ SPS5D = 1/2; % 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
+ 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= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
+ %% 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 = 0.2; %
+ 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;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+ if ((RERUN || isempty(data.NU_EVOL)) && RUN)
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+
+ % Load results after trying to run
+ 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);
+% g_ky(i,j,:) = g_ky;
+%
+% g_avg(i,j,:) = lg.avg_g;
+% g_std(i,j,:) = lg.std_g;
+
+ [~,it1] = min(abs(data.Ts3D-0.5*data.Ts3D(end))); % start of the measurement time window
+ [~,it2] = min(abs(data.Ts3D-1.0*data.Ts3D(end))); % end of ...
+ field = 0;
+ field_t = 0;
+ for ik = 1:NY/2+1
+ field = squeeze(sum(abs(data.PHI),3)); % take the sum over z
+ field_t = squeeze(field(ik,1,:)); % take the kx =0, ky = ky mode only
+ to_measure = log(field_t);
+ gr = polyfit(data.Ts3D(it1:it2),to_measure(it1:it2),1);
+ g_ky(i,j,ik) = gr(1);
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data.ky(ikmax)); disp(msg);
+
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 1
+%% Study of the peak growth rate
+figure
+
+y_ = g_ky;
+e_ = 0.05;
+
+% filter to noisy data
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+
+[y_,idx_] = max(g_ky,[],3);
+for i = 1:numel(idx_)
+ e_ = g_std(:,:,idx_(i));
+end
+
+colors_ = lines(numel(NU_a));
+subplot(121)
+for i = 1:numel(NU_a)
+% errorbar(P_a,y_(i,:),e_(i,:),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
+% 'color',colors_(i,:));
+ plot(P_a,y_(i,:),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['$\nu=$',num2str(NU_a(i))],...
+ 'color',colors_(i,:));
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
+legend('show'); xlabel('$P$, $J=P/2$'); ylabel('$\gamma$');
+
+colors_ = jet(numel(P_a));
+subplot(122)
+for j = 1:numel(P_a)
+% errorbar(NU_a,y_(:,j),e_(:,j),...
+% 'LineWidth',1.2,...
+% 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+% 'color',colors_(j,:));
+ plot(NU_a,y_(:,j),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',['(',num2str(P_a(j)),',',num2str(J_a(j)),')'],...
+ 'color',colors_(j,:));
+ hold on;
+end
+title(['$\kappa_T=$',num2str(K_Ti),' $k_y=k_y^{max}$']);
+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/analysis_gene.m b/wk/analysis_gene.m
index 5e609cbb8b59a9a1015b6aa4514de58665221467..42da732b67b4313c71984d2ed417b3a6936c2499 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -13,7 +13,7 @@ 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_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/';
@@ -27,7 +27,8 @@ folder = '/misc/gene_results/Z-pinch/HP_fig_2a_gyroLES/';
% folder = '/misc/gene_results/Z-pinch/kN_2.5_HD_transport_spectrum_01/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/196x96x20x32x16_02/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
+% folder = '/misc/gene_results/CBC/256x96x24x32x12/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
@@ -41,6 +42,9 @@ folder = '/misc/gene_results/Z-pinch/HP_fig_2a_gyroLES/';
% folder = '/misc/gene_results/CBC/KT_5.3_192x96x24x30x16_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
+% folder = '/misc/gene_results/linear_miller_CBC/hdf5_miller_s0_adiabe/';
+folder = '/misc/gene_results/linear_miller_CBC/hdf5_salpha_s0_adiabe/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -70,15 +74,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 = 'xy';
+options.PLAN = 'kxz';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [1:10];
+options.TIME = [40 55 70];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -89,11 +93,11 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'v_y';
% options.NAME = '\Gamma_x';
-options.NAME = 'n_i';
-options.PLAN = 'xy';
+% options.NAME = 'n_i';
+options.PLAN = 'kyz';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -122,7 +126,7 @@ subplot(312)
xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
subplot(313)
- for i = 11:16
+ for i = [11 12 14]
plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
@@ -131,9 +135,12 @@ end
if 0
%% Show f_i(vpar,mu)
-options.times = 100:200;
+options.times = 70;
options.specie = 'i';
-options.PLT_FCT = 'contour';
+% options.PLT_FCT = 'contour';
+% options.PLT_FCT = 'contourf';
+% options.PLT_FCT = 'surf';
+options.PLT_FCT = 'surfvv';
options.folder = folder;
options.iz = 'avg';
options.FIELD = '<f_>';
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index f787c1c23b19af625bc453902c0ae45d67b87358..50b84ffac55fc35e80ebaa84e88bd5f6f36ebd87 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -69,13 +69,13 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'yz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [6000:1:9000];
+options.TIME = [1:1:9000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -85,9 +85,9 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
-options.AXISEQUAL = 1;
+options.AXISEQUAL = 0;
options.NAME = '\phi';
% options.NAME = '\psi';
% options.NAME = 'n_i';
@@ -95,11 +95,11 @@ options.NAME = '\phi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
-options.NAME = 'f_i';
-options.PLAN = 'sx';
+options.PLAN = 'kxz';
+% options.NAME = 'f_i';
+% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [5 20 50 100 150];
+options.TIME = [100 250 300];
options.RESOLUTION = 256;
data.a = data.EPS * 2e3;
@@ -124,12 +124,15 @@ 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);
+options.XPERP = linspace( 0,sqrt(6),16).^2;
+% options.SPAR = gene_data.vp';
+% options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [100 150];
-options.PLT_FCT = 'contour';
+options.T = [100];
+% options.PLT_FCT = 'contour';
+% options.PLT_FCT = 'contourf';
+options.PLT_FCT = 'surfvv';
options.ONED = 0;
options.non_adiab = 0;
options.SPECIE = 'i';
@@ -142,11 +145,11 @@ if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
options.P2J = 0;
-options.ST = 1;
+options.ST = 0;
options.PLOT_TYPE = 'space-time';
options.NORMALIZED = 0;
options.JOBNUM = 0;
-options.TIME = [1000];
+options.TIME = [200:500];
options.specie = 'i';
options.compz = 'avg';
fig = show_moments_spectrum(data,options);
@@ -212,8 +215,8 @@ end
if 0
%% Metric infos
-options.SHOW_FLUXSURF = 1;
-options.SHOW_METRICS = 0;
+options.SHOW_FLUXSURF = 0;
+options.SHOW_METRICS = 1;
fig = plot_metric(data,options);
end
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 4dc2a135404afd4df84ddf45b6f0ecf675c69de6..34e192e05ff96dfc1173546594a281dc2db60026 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -31,7 +31,7 @@ gyacomodir = '/home/ahoffman/gyacomo/';
%% CBC
% resdir = 'CBC/64x32x16x5x3';
% resdir = 'CBC/64x128x16x5x3';
-% resdir = 'CBC/old/128x64x16x5x3';
+resdir = 'CBC/old/128x64x16x5x3';
% resdir = 'CBC/96x96x16x3x2_Nexc_6';
% resdir = 'CBC/128x96x16x3x2_Nexc_0';
% resdir = 'CBC/old/192x96x24x13x7';
@@ -45,7 +45,7 @@ gyacomodir = '/home/ahoffman/gyacomo/';
% 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/64x32x16x3x2_Nexc_0_change_dt';
% resdir = 'CBC/128x96x16x3x2_Nexc_0';
% resdir = 'CBC/128x96x16x3x2_Nexc_0';
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 039b3f4950f627221779d6ec635b41b755eb296d..a54a321ae08163ee647bc6e9763c76f47c43cdab 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -3,15 +3,19 @@
% 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
+
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0';
SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
-addpath(genpath('../matlab')) % ... add
default_plots_options
-gyacomodir = '/home/ahoffman/gyacomo/';
-% EXECNAME = 'gyacomo_1.0';
% EXECNAME = 'gyacomo_dbg';
-EXECNAME = 'gyacomo';
+EXECNAME = 'gyacomo_alpha';
+% EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
@@ -19,26 +23,26 @@ 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.22; % ele Density '''
-K_Te = 6.96; % ele Temperature '''
-K_Ni = 2.22; % ion Density gradient drive
+K_Ne = 1*2.22; % ele Density '''
+K_Te = 1*6.96; % ele Temperature '''
+K_Ni = 1*2.22; % ion Density gradient drive
K_Ti = 6.96; % 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 = 0; % 1: kinetic electrons, 2: adiabatic electrons
-BETA = 0.0001; % electron plasma beta
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
%% GRID PARAMETERS
-P = 4;
+P = 8;
J = P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 2; % real space x-gridpoints
-NY = 12; % '' y-gridpoints
+NX = 8; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.8; % 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)
+LY = 2*pi/0.4; % Size of the squared frequency domain
+NZ = 24; % 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))
@@ -47,19 +51,21 @@ GEOMETRY= 's-alpha';
% GEOMETRY= 'miller';
EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
-SHEAR = 0.0; % magnetic shear
+SHEAR = 0.8; % magnetic shear
KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
-PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+% PARALLEL_BC = 'dirichlet'; %'dirichlet','periodic','shearless','disconnected'
+PARALLEL_BC = 'periodic'; %'dirichlet','periodic','shearless','disconnected'
SHIFT_Y = 0.0;
%% TIME PARMETERS
-TMAX = 100; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 50; % Maximal time unit
+% DT = 1e-2; % Time step
+DT = 1e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = -1; % Sampling per time unit for 2D arrays
-SPS3D = 2; % Sampling per time unit for 2D arrays
-SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPS3D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/2; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
JOB2LOAD= -1;
%% OPTIONS
@@ -73,7 +79,7 @@ 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
+INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
NUMERICAL_SCHEME = 'RK4'; % RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5
%% OUTPUTS
W_DOUBLE = 1;
@@ -90,7 +96,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 1.0; %
+MU_Z = 0.2; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
@@ -106,19 +112,22 @@ setup
% Run linear simulation
if RUN
% 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 4 ',gyacomodir,'bin/',EXECNAME,' 2 2 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 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
end
%% Load results
%%
filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [gyacomodir,'results/',filename,'/'];
+LOCALDIR = [gyacomodir,'results/',filename,'/'];
+FIGDIR = LOCALDIR;
% Load outputs from jobnummin up to jobnummax
JOBNUMMIN = 00; JOBNUMMAX = 01;
data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+data.FIGDIR = ['../results/',SIMID,'/',PARAMS,'/'];
%% Short analysis
if 0
@@ -133,7 +142,7 @@ 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
+if 0
%% Ballooning plot
options.time_2_plot = [10 50];
options.kymodes = [0.1 0.2 0.4];
@@ -148,9 +157,10 @@ if 0
options.P2J = 1;
options.ST = 1;
options.PLOT_TYPE = 'space-time';
-% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-1D';ls
+
% options.PLOT_TYPE = 'Tavg-2D';
-options.NORMALIZED = 0;
+options.NORMALIZED = 1;
options.JOBNUM = 0;
options.TIME = [0 50];
options.specie = 'i';
@@ -175,15 +185,17 @@ if 1
options.NORMALIZED = 0;
options.K2PLOT = 1;
options.TIME = [0:1000];
+% options.TIME = [0.5:0.01:1].*data.Ts3D(end);
options.NMA = 1;
-options.NMODES = 16;
+options.NMODES = 32;
options.iz = 'avg';
+options.ik = 1;
fig = mode_growth_meter(data,options);
save_figure(data,fig,'.png')
end
-if 1
+if 0
%% RH TEST
ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));