diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
index 6bb81548c6007ccde0c3991a2101434c12b7ea95..6c0afb15853c1cde9e9410f546d2b905709035a6 100644
--- a/matlab/assemble_cosolver_matrices.m
+++ b/matlab/assemble_cosolver_matrices.m
@@ -1,6 +1,6 @@
codir = '/home/ahoffman/cosolver/';
-matname= 'gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5';
-matdir = [codir,'SGGK_tau1e-2_P4_J2_dk_0.05_kpm_5.0/matrices/'];
+matname= 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8';
+matdir = [codir,'LDGKii_P16_J9_NFLR_8/matrices/'];
kperp = load([matdir,'kperp.in']);
Nkp = numel(kperp);
outdir = '/home/ahoffman/gyacomo/iCa/';
@@ -13,34 +13,34 @@ end
Nmax = numel(kperp);
for n_=1:Nmax
n_string = sprintf('%5.5d',n_); disp(n_string);
- filename = ['ei.',n_string,'.h5'];
+ filename = ['self.',n_string,'.h5'];
if(exist([matdir,filename])>0)
% Load matrices and write them in one h5 file
- olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
- newdname = ['/',sprintf('%5.5d',n_-1),olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ceipj/CeipjT'; infile = [matdir,'ei.',n_string,'.h5'];
- newdname = ['/',sprintf('%5.5d',n_-1),olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ciepj/CiepjT'; infile = [matdir,'ie.',n_string,'.h5'];
- newdname = ['/',sprintf('%5.5d',n_-1),olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Ciepj/CiepjF'; infile = [matdir,'ie.',n_string,'.h5'];
- newdname = ['/',sprintf('%5.5d',n_-1),olddname];
- load_write_mat(infile,olddname,outfile,newdname);
-
- olddname = '/Caapj/Ceepj'; infile = [matdir,'self.',n_string,'.h5'];
- newdname = ['/',sprintf('%5.5d',n_-1),olddname];
- mat_ee = load_write_mat(infile,olddname,outfile,newdname);
+ % olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
+ % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ % load_write_mat(infile,olddname,outfile,newdname);
+ %
+ % olddname = '/Ceipj/CeipjT'; infile = [matdir,'ei.',n_string,'.h5'];
+ % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ % load_write_mat(infile,olddname,outfile,newdname);
+ %
+ % olddname = '/Ciepj/CiepjT'; infile = [matdir,'ie.',n_string,'.h5'];
+ % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ % load_write_mat(infile,olddname,outfile,newdname);
+ %
+ % olddname = '/Ciepj/CiepjF'; infile = [matdir,'ie.',n_string,'.h5'];
+ % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ % load_write_mat(infile,olddname,outfile,newdname);
+ %
+ % olddname = '/Caapj/Ceepj'; infile = [matdir,'self.',n_string,'.h5'];
+ % newdname = ['/',sprintf('%5.5d',n_-1),olddname];
+ % mat_ee = load_write_mat(infile,olddname,outfile,newdname);
olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
newdname = ['/',sprintf('%5.5d',n_-1),olddname];
mat_ii = load_write_mat(infile,olddname,outfile,newdname);
-
+ mat_ee = 0;
% to verify symmetry
verif = max(abs(imag(eig(mat_ee)))) + max(abs(imag(eig(mat_ii))));
if(verif>0) disp(['Warning, non sym matrix at ', n_string]); end;
@@ -51,16 +51,16 @@ for n_=1:Nmax
break
end
end
-olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',n_string,'.h5'];
-Pmaxe = h5readatt(infile,olddname,'Pmaxe');
-Jmaxe = h5readatt(infile,olddname,'Jmaxe');
+olddname = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
+% Pmaxe = h5readatt(infile,olddname,'Pmaxe');
+% Jmaxe = h5readatt(infile,olddname,'Jmaxe');
Pmaxi = h5readatt(infile,olddname,'Pmaxi');
Jmaxi = h5readatt(infile,olddname,'Jmaxi');
-dims_e = [0 0];
-dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
+% dims_e = [0 0];
+% dims_e(1)= Pmaxe; dims_e(2) = Jmaxe;
dims_i = [Pmaxi; Jmaxi];
-h5create(outfile,'/dims_e',numel(dims_e));
-h5write (outfile,'/dims_e',dims_e);
+% h5create(outfile,'/dims_e',numel(dims_e));
+% h5write (outfile,'/dims_e',dims_e);
h5create(outfile,'/dims_i',numel(dims_i));
h5write (outfile,'/dims_i',dims_i);
%
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 73f72a925f221167857f3fabf6ab15a1e5f3312b..191cdf1a6f2650d4a6cc4f5c5437b4f7d9e298d4 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -1,7 +1,6 @@
function [DATA] = load_params(DATA,filename)
DATA.inputs.CO = h5readatt(filename,'/data/input/coll','CO');
- DATA.inputs.K_N = h5readatt(filename,'/data/input/ions','k_N');
- DATA.inputs.K_T = h5readatt(filename,'/data/input/ions','k_T');
+
DATA.inputs.Q0 = h5readatt(filename,'/data/input/geometry','q0');
DATA.inputs.EPS = h5readatt(filename,'/data/input/geometry','eps');
DATA.inputs.SHEAR = h5readatt(filename,'/data/input/geometry','shear');
@@ -38,10 +37,12 @@ function [DATA] = load_params(DATA,filename)
DATA.inputs.MUz = h5readatt(filename,'/data/input/model','mu_z');
DATA.inputs.LINEARITY = h5readatt(filename,'/data/input/model','LINEARITY');
DATA.inputs.BETA = h5readatt(filename,'/data/input/model','beta');
- DATA.inputs.TAU_E = h5readatt(filename,'/data/input/model','tau_e');
+ DATA.inputs.TAU_I = h5readatt(filename,'/data/input/model','tau_i');
DATA.inputs.HYP_V = h5readatt(filename,'/data/input/model','HYP_V');
DATA.inputs.K_cB = h5readatt(filename,'/data/input/model','k_cB');
DATA.inputs.K_gB = h5readatt(filename,'/data/input/model','k_gB');
+ DATA.inputs.ADIAB_E = h5readatt(filename,'/data/input/model','ADIAB_E') == 'y';
+ DATA.inputs.ADIAB_I = h5readatt(filename,'/data/input/model','ADIAB_I') == 'y';
DATA.inputs.W_GAMMA = h5readatt(filename,'/data/input/diag_par','write_gamma') == 'y';
DATA.inputs.W_PHI = h5readatt(filename,'/data/input/diag_par','write_phi') == 'y';
@@ -56,6 +57,13 @@ function [DATA] = load_params(DATA,filename)
DATA.inputs.K_N = zeros(1,DATA.inputs.Na);
DATA.inputs.K_T = zeros(1,DATA.inputs.Na);
spnames = {'ions','electrons'};
+ if(DATA.inputs.ADIAB_E)
+ spnames = {spnames{1}};
+ end
+ if(DATA.inputs.ADIAB_I)
+ spnames = {spnames{2}};
+ end
+
for ia=1:DATA.inputs.Na
spdata = ['/data/input/',spnames{ia}];
DATA.inputs.sigma(ia) = h5readatt(filename,spdata,'sigma');
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index a37b3b714f623306f9193552b446ca552916f681..0b8eaea271fb1efa6af82169e6a62a7dc47ddab5 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -2,20 +2,19 @@ function [FIG] = plot_ballooning(data,options)
FIG.fig = figure; iplot = 1;
[~,it0] = min(abs(data.Ts3D - options.time_2_plot(1)));
[~,it1] = min(abs(data.Ts3D - options.time_2_plot(end)));
- [~,ikyarray] = min(abs(data.ky - options.kymodes));
-% phi_real=mean(real(data.PHI(:,:,:,it0:it1)),4);
-% phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
+ [~,ikyarray] = min(abs(data.grids.ky - options.kymodes));
+ ikyarray = unique(ikyarray);
phi_real=real(data.PHI(:,:,:,it1));
phi_imag=imag(data.PHI(:,:,:,it1));
ncol = 1;
- if data.BETA > 0
+ if data.inputs.BETA > 0
psi_real=real(data.PSI(:,:,:,it1));
psi_imag=imag(data.PSI(:,:,:,it1));
ncol = 2;
end
% Apply ballooning transform
- if(data.Nkx > 1)
- nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2));
+ if(data.grids.Nkx > 1)
+ nexc = round(data.grids.ky(2)*data.inputs.SHEAR*2*pi/data.grids.kx(2));
else
nexc = 1;
end
@@ -34,7 +33,7 @@ function [FIG] = plot_ballooning(data,options)
ordered_ikx = [(Np_+1):Nkx 1:Np_];
end
try
- idx=data.kx./data.kx(2);
+ idx=data.grids.kx./data.grids.kx(2);
catch
idx=0;
end
@@ -54,11 +53,11 @@ function [FIG] = plot_ballooning(data,options)
end
% Define ballooning angle
- coordz = data.z;
+ coordz = data.grids.z;
for i_ =1: Nkx
for iz=1:dims(3)
ii = dims(3)*(i_-1) + iz;
- b_angle(ii) =coordz(iz) + 2*pi*idx(i_)./is;
+ b_angle(ii) =coordz(iz) + 2*pi*data.grids.Npol*idx(i_)./is;
end
end
@@ -73,6 +72,7 @@ function [FIG] = plot_ballooning(data,options)
phib_norm = phib / normalization;
phib_real_norm = real( phib_norm);
phib_imag_norm = imag( phib_norm);
+ %
subplot(numel(ikyarray),ncol,ncol*(iplot-1)+1)
plot(b_angle/pi, phib_real_norm,'-b'); hold on;
plot(b_angle/pi, phib_imag_norm ,'-r');
@@ -80,18 +80,20 @@ function [FIG] = plot_ballooning(data,options)
legend('real','imag','norm')
xlabel('$\chi / \pi$')
ylabel('$\phi_B (\chi)$');
- title(['$k_y=',sprintf('%2.2f',data.ky(iky)),...
+ title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
-
- if data.BETA > 0
+ for i_ = 2*[1:data.grids.Nkx]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ end
+ if data.inputs.BETA > 0
psib_real = zeros( Nkx*dims(3) ,1);
psib_imag = psib_real;
for i_ =1:Nkx
- start_ = (i_ -1)*dims(3) +1;
+ start_ = (i_-1)*dims(3) +1;
end_ = i_*dims(3);
ikx = ordered_ikx(i_);
- psib_real(start_:end_) = psi_real(iky,ikx,:);
- psib_imag(start_:end_) = psi_imag(iky,ikx,:);
+ psib_real(start_:end_) = psi_real(iky,ikx,:);
+ psib_imag(start_:end_) = psi_imag(iky,ikx,:);
end
psib = psib_real(:) + 1i * psib_imag(:);
psib_norm = psib / normalization;
@@ -102,9 +104,12 @@ function [FIG] = plot_ballooning(data,options)
plot(b_angle/pi, psib_imag_norm ,'-r');
plot(b_angle/pi, abs(psib_norm),'-k');
legend('real','imag','norm')
+ for i_ = 2*[1:data.grids.Nkx]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ end
xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$');
- title(['$k_y=',sprintf('%2.2f',data.ky(iky)),...
+ title(['$k_y=',sprintf('%2.2f',data.grids.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
end
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 8636462485f2973228490eee455185cbb8a4ae60..8a813cff4501b6cd5cff46fccc8874a82b8e7e0a 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -108,9 +108,10 @@ mfns = {...
'/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';...
'/home/ahoffman/gyacomo/iCa/gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';...
'/home/ahoffman/gyacomo/iCa/gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';...
- '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';...
- '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
- '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';...
+ % '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P11_J7_dk_5e-2_km_5.0_NFLR_12.h5';...
+ % '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-3_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
+ % '/home/ahoffman/gyacomo/iCa/gk_sugama_tau1e-2_P4_J2_dk_5e-2_km_5.0_NFLR_5.h5';...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5';...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_4_J_2_N_75_kpm_5.0.h5';...
};
@@ -125,21 +126,13 @@ CONAME_A = {...
'FC 4 2 NFLR 12'; ...
'LD 10 5 NFLR 12'; ...
'LD 11 7 NFLR 16'; ...
- 'SG 11 7 NFLR 12, tau 1e-3'; ...
- 'SG 4 2 NFLR 5, tau 1e-3'; ...
- 'SG 4 2 NFLR 5, tau 1e-2'; ...
+ 'LDii 16 9 NFLR 8'; ...
+ % 'SG 11 7 NFLR 12, tau 1e-3'; ...
+ % 'SG 4 2 NFLR 5, tau 1e-3'; ...
+ % 'SG 4 2 NFLR 5, tau 1e-2'; ...
% 'Hacked SG A';...
% 'Hacked SG B';...
};
-TAU_A = [1;...
- 1;...
- 1;...
- 1;...
- 1;...
- 1e-3;...
- 1e-3;...
- 1e-2;...
- ];
figure
for j_ = 1:numel(mfns)
mat_file_name = mfns{j_}; disp(mat_file_name);
@@ -154,9 +147,9 @@ for j_ = 1:numel(mfns)
wmax(idx_+1) = max(abs(imag(eig(MAT))));
end
subplot(121)
- plot(kp_a,gmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on;
+ plot(kp_a,gmax,'DisplayName',CONAME_A{j_}); hold on;
subplot(122)
- plot(kp_a,wmax*TAU_A(j_),'DisplayName',CONAME_A{j_}); hold on;
+ plot(kp_a,wmax,'DisplayName',CONAME_A{j_}); hold on;
end
subplot(121)
legend('show'); grid on;
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index f03fa6be93901edc99a98408b226eaa639275bf8..d1caeced8f68f96b0da031403ae1ce41fe642f3d 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -5,42 +5,43 @@ function [ fig, arrays ] = plot_metric( data, options )
names = {'gxx','gxy','gxz','gyy','gyz','gzz',...
'hatB', 'dBdx', 'dBdy', 'dBdz',...
'Jacobian','hatR','hatZ','gradz_coeff'};
-geo_arrays = zeros(2,data.Nz,numel(names));
+geo_arrays = zeros(data.grids.Nz,numel(names));
for i_ = 1:numel(names)
namae = names{i_};
- geo_arrays(:,:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
+ geo_arrays(:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end
NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
if NPLOT > 0
fig = figure;
if options.SHOW_METRICS
- subplot(311)
+ subplot(3,NPLOT,1*NPLOT)
for i = 1:6
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
+ xlim([min(data.grids.z),max(data.grids.z)]); legend('show'); title('GYACOMO geometry');
- subplot(312)
+ subplot(3,NPLOT,2*NPLOT)
for i = 7:10
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(data.z),max(data.z)]); legend('show');
+ xlim([min(data.grids.z),max(data.grids.z)]); legend('show');
- subplot(313)
+ subplot(3,NPLOT,3*NPLOT)
for i = 11:14
- plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
+ plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(data.z),max(data.z)]); legend('show');
+ xlim([min(data.grids.z),max(data.grids.z)]); legend('show');
end
if options.SHOW_FLUXSURF
- R = [geo_arrays(1,:,12),geo_arrays(1,1,12)];
- Z = [geo_arrays(1,:,13),geo_arrays(1,1,13)];
+ subplot(1,2,1)
+ R = [geo_arrays(:,12);geo_arrays(1,12)];
+ Z = [geo_arrays(:,13);geo_arrays(1,13)];
plot(R,Z);
axis equal
end
end
%outputs
-arrays = squeeze(geo_arrays(1,:,:));
+arrays = squeeze(geo_arrays(:,:));
end
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index dca0292767ff0bc27b4e710544b2db1222c6d708..01bc00d7c5bf341ce6f262b9d44bd380250eeb03 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -1,18 +1,18 @@
function [ FIGURE ] = show_geometry(DATA,OPTIONS)
% filtering Z pinch and torus
-if DATA.Nz > 1 % Torus flux tube geometry
- Nturns = floor(abs(DATA.z(1)/pi));
- Nptor = ceil(DATA.Nz*2/Nturns); Tgeom = 1;
+if DATA.grids.Nz > 1 % Torus flux tube geometry
+ Nturns = floor(abs(DATA.grids.z(1)/pi));
+ Nptor = ceil(DATA.grids.Nz*2/Nturns); Tgeom = 1;
a = DATA.EPS; % Torus minor radius
R = 1.; % Torus major radius
- q = (DATA.Nz>1)*DATA.Q0; % Flux tube safety factor
+ q = (DATA.grids.Nz>1)*DATA.inputs.Q0; % Flux tube safety factor
theta = linspace(-pi, pi, Nptor) ; % Poloidal angle
phi = linspace(0., 2.*pi, Nptor) ; % Toroidal angle
[t, p] = meshgrid(phi, theta);
x_tor = (R + a.*cos(p)) .* cos(t);
y_tor = (R + a.*cos(p)) .* sin(t);
z_tor = a.*sin(p);
- DIMENSIONS = [50 50 1200 600];
+ DIMENSIONS = [600 600 1200 600];
else % Zpinch geometry
Nturns = 0.1; Tgeom = 0;
a = 0.7; % Torus minor radius
@@ -48,7 +48,7 @@ yZ = @(z) bX(z).*xY(z)-bY(z).*xX(z)./sqrt((bY(z).*xZ(z)-bZ(z).*xY(z)).^2+(bZ(z)
yvec= @(z) [yX(z); yY(z); yZ(z)];
% Plot high res field line
% Planes plot
-theta = DATA.z ; % Poloidal angle
+theta = DATA.grids.z ; % Poloidal angle
% Plot x,y,bvec at these points
scale = 0.10;
@@ -56,11 +56,11 @@ scale = 0.10;
OPTIONS.POLARPLOT = 0;
OPTIONS.PLAN = 'xy';
r_o_R = DATA.rho_o_R;
-[X,Y] = meshgrid(r_o_R*DATA.x,r_o_R*DATA.y);
+[X,Y] = meshgrid(r_o_R*DATA.grids.x,r_o_R*DATA.grids.y);
max_ = 0;
-FIELDS = zeros(DATA.Ny,DATA.Nx,DATA.Nz);
+FIELDS = zeros(DATA.grids.Ny,DATA.grids.Nx,DATA.grids.Nz);
-FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.PARAMS]; set(gcf, 'Position', DIMENSIONS)
+FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.params_string]; set(gcf, 'Position', DIMENSIONS)
for it_ = 1:numel(OPTIONS.TIME)
subplot(1,numel(OPTIONS.TIME),it_)
%plot magnetic geometry
@@ -86,7 +86,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
end
end
%plot vector basis
- theta = DATA.z ; % Poloidal angle
+ theta = DATA.grids.z ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*xX(theta),scale*xY(theta),scale*xZ(theta),0,'r');
quiver3(Xfl(theta),Yfl(theta),Zfl(theta),scale*yX(theta),scale*yY(theta),scale*yZ(theta),0,'g');
@@ -101,7 +101,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
if (tmp > max_) max_ = tmp;
end
for iz = OPTIONS.PLANES
- z_ = DATA.z(iz);
+ z_ = DATA.grids.z(iz);
Xp = Xfl(z_) + X*xX(z_) + Y*yX(z_);
Yp = Yfl(z_) + X*xY(z_) + Y*yY(z_);
Zp = Zfl(z_) + X*xZ(z_) + Y*yZ(z_);
@@ -110,7 +110,7 @@ subplot(1,numel(OPTIONS.TIME),it_)
colormap(bluewhitered);
% caxis([-1,1]);
end
- if DATA.Nz == 1
+ if DATA.grids.Nz == 1
xlim([0.65 0.75]);
ylim([-0.1 0.1]);
zlim([-0.2 0.2]);
diff --git a/matlab/setup.m b/matlab/setup.m
index 9e0490c8a1902dee8396db7962a8a94a9c9b88a7..a5b6f420fb80ca48b15bc693a96a61445b9beb9f 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -15,17 +15,28 @@ GEOM.geom = ['''',GEOMETRY,''''];
GEOM.q0 = Q0; % q factor
GEOM.shear = SHEAR; % shear
GEOM.eps = EPS; % inverse aspect ratio
-GEOM.kappa = KAPPA; % elongation
-GEOM.delta = DELTA; % triangularity
-GEOM.zeta = ZETA; % squareness
+GEOM.kappa = KAPPA; % elongation
+GEOM.s_kappa = S_KAPPA;
+GEOM.delta = DELTA; % triangularity
+GEOM.s_delta = S_DELTA;
+GEOM.zeta = ZETA; % squareness
+GEOM.s_zeta = S_ZETA;
GEOM.parallel_bc = ['''',PARALLEL_BC,''''];
GEOM.shift_y = SHIFT_Y;
GEOM.Npol = NPOL;
+if PB_PHASE; GEOM.PB_PHASE = '.true.'; else; GEOM.PB_PHASE = '.false.';end;
% Model parameters
MODEL.LINEARITY = ['''',LINEARITY,''''];
+try
+ RM_LD_T_EQ;
+catch
+ RM_LD_T_EQ = 0;
+end
if RM_LD_T_EQ; MODEL.RM_LD_T_EQ = '.true.'; else; MODEL.RM_LD_T_EQ = '.false.'; end;
MODEL.Na = NA;
if ADIAB_E; MODEL.ADIAB_E = '.true.'; else; MODEL.ADIAB_E = '.false.';end;
+if ADIAB_I; MODEL.ADIAB_I = '.true.'; else; MODEL.ADIAB_I = '.false.';end;
+if MHD_PD; MODEL.MHD_PD = '.true.'; else; MODEL.MHD_PD = '.false.';end;
MODEL.beta = BETA;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
@@ -44,7 +55,6 @@ MODEL.sigma_i = 1.0;
% charge q_a/e
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
-if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
MODEL.K_Ni = K_Ni;
MODEL.K_Ne = K_Ne;
@@ -71,7 +81,8 @@ switch CO
case 'LR'
COLL.mat_file = 'gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
case 'LD'
- COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
+ % COLL.mat_file = 'gk_landau_P10_J5_dk_5e-2_km_2.0_NFLR_12.h5';
+ COLL.mat_file = 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8.h5';
% COLL.mat_file = 'gk_landau_P11_J7_dk_5e-2_km_2.0_NFLR_16.h5';
% COLL.mat_file = 'gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
% COLL.mat_file = 'LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 67520b5a0a8516b5e631ce9a0cf5a4623fb7f1ad..555931e62f8ddb7193996542ce716d03ac91c811 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -29,11 +29,15 @@ fprintf(fid,[' q0 = ', num2str(GEOM.q0),'\n']);
fprintf(fid,[' shear = ', num2str(GEOM.shear),'\n']);
fprintf(fid,[' eps = ', num2str(GEOM.eps),'\n']);
fprintf(fid,[' kappa = ', num2str(GEOM.kappa),'\n']);
+fprintf(fid,['s_kappa = ', num2str(GEOM.s_kappa),'\n']);
fprintf(fid,[' delta = ', num2str(GEOM.delta),'\n']);
+fprintf(fid,['s_delta = ', num2str(GEOM.s_delta),'\n']);
fprintf(fid,[' zeta = ', num2str(GEOM.zeta),'\n']);
+fprintf(fid,['s_zeta = ', num2str(GEOM.s_zeta),'\n']);
fprintf(fid,[' parallel_bc = ', GEOM.parallel_bc,'\n']);
fprintf(fid,[' shift_y = ', num2str(GEOM.shift_y),'\n']);
fprintf(fid,[' Npol = ', num2str(GEOM.Npol),'\n']);
+fprintf(fid,[' PB_PHASE= ', GEOM.PB_PHASE,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&OUTPUT_PAR\n');
@@ -69,6 +73,9 @@ fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
fprintf(fid,[' ADIAB_E = ', MODEL.ADIAB_E,'\n']);
+fprintf(fid,[' ADIAB_I = ', MODEL.ADIAB_I,'\n']);
+fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
+fprintf(fid,[' MHD_PD = ', MODEL.MHD_PD,'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&CLOSURE_PAR\n');
@@ -78,18 +85,19 @@ fprintf(fid,[' nonlinear_closure=',CLOSURE.nonlinear_closure,'\n']);
fprintf(fid,[' nmax =',num2str(CLOSURE.nmax),'\n']);
fprintf(fid,'/\n');
-fprintf(fid,'&SPECIES\n');
-fprintf(fid, ' name_ = ions \n');
-fprintf(fid,[' tau_ = ', num2str(MODEL.tau_i),'\n']);
-fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_i),'\n']);
-fprintf(fid,[' q_ = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ni),'\n']);
-fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Ti),'\n']);
-fprintf(fid,'/\n');
-
-if(MODEL.Na > 1)
- fprintf(fid,'&SPECIES\n');
- fprintf(fid, ' name_ = electrons');
+if(strcmp(MODEL.ADIAB_I,'.false.'))
+ fprintf(fid,'&SPECIES\n');
+ fprintf(fid, ' name_ = ions \n');
+ fprintf(fid,[' tau_ = ', num2str(MODEL.tau_i),'\n']);
+ fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_i),'\n']);
+ fprintf(fid,[' q_ = ', num2str(MODEL.q_i),'\n']);
+ fprintf(fid,[' K_N_ = ', num2str(MODEL.K_Ni),'\n']);
+ fprintf(fid,[' K_T_ = ', num2str(MODEL.K_Ti),'\n']);
+ fprintf(fid,'/\n');
+end
+if(strcmp(MODEL.ADIAB_E,'.false.'))
+ fprintf(fid,'&SPECIES\n');
+ fprintf(fid, ' name_ = electrons \n');
fprintf(fid,[' tau_ = ', num2str(MODEL.tau_e),'\n']);
fprintf(fid,[' sigma_ = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' q_ = ', num2str(MODEL.q_e),'\n']);
diff --git a/wk/load_metadata_scan.m b/wk/CBC_load_metadata_scan.m
similarity index 87%
rename from wk/load_metadata_scan.m
rename to wk/CBC_load_metadata_scan.m
index 040038705b053b23e0832d9e2b4a01c5c9b5ef8b..6c9e87d230c0931b580c051044d34ea7c8d58b75 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/CBC_load_metadata_scan.m
@@ -60,15 +60,26 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'p2_linear_new/8x24_ky_0.3_P_16_J_8_kT_2.5_4.5_nu_0.001_0.1_DGGK.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_5.3.mat';
%% ky pj scan
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_64_DGDK_0.001_kT_6.96.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_48_DGDK_0.001_kT_6.96.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_DGGK_0.01_kT_5.3.mat';
% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_16_SGGK_0.01_kT_5.3.mat';
-datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_10_LDGK_0.01_kT_5.3.mat';
-%% Data for the paper :
-% Collisionless kT threshold
-% datafname = 'p2_linear_new/8x24_ky_0.3_kT_3_6.96_P_2_30_DGDK_0.001.mat';
+% datafname = 'p2_linear_new/8x24_ky_0.05_1_P_2_10_LDGK_0.01_kT_5.3.mat';
+
+% datafname = 'p2_linear_new/8x24_P_4_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_DGGK_nu_0.005_0.05_kT_5.3.mat';
+
+% datafname = 'p2_linear_new/8x24_P_4_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
+datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_SGGK_nu_0.005_0.05_kT_5.3.mat';
+
+% datafname = 'p2_linear_new/8x24_P_4_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_8_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_10_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
+% datafname = 'p2_linear_new/8x24_P_16_ky_0.05_1_LDGK_nu_0.005_0.05_kT_5.3.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
@@ -104,14 +115,15 @@ if 1
figure
colors_ = jet(numel(d.s2));
for i = 1:numel(d.s2)
-% plot(d.s1,d.data(:,i),'s-',...
-% 'LineWidth',2.0,...
-% 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
-% 'color',colors_(i,:));
- errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
- 'LineWidth',2.0,...
- 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
- 'color',colors_(i,:));
+ % plot(d.s1,d.data(:,i),'s-',...
+ plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ 'color',colors_(i,:));
+ % errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
+ % 'LineWidth',2.0,...
+ % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ % 'color',colors_(i,:));
hold on;
end
xlabel(d.s1name); ylabel(d.dname);title(d.title);
diff --git a/wk/ETG_adiab_i.m b/wk/ETG_adiab_i.m
new file mode 100644
index 0000000000000000000000000000000000000000..ce0bdb1087fb9200b890b212b8ccce7d468d54fb
--- /dev/null
+++ b/wk/ETG_adiab_i.m
@@ -0,0 +1,162 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'ETG_adiab_i'; % Name of the simulation
+RUN = 1; % To run or just to load
+default_plots_options
+% EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
+EXECNAME = 'gyacomo23_debug'; % single precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.0; % Collision frequency
+TAU = 1; % e/i temperature ratio
+K_Ne = 0; % ele Density '''
+K_Te = 10; % ele Temperature '''
+K_Ni = 0; % ion Density gradient drive
+K_Ti = 0; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 1; % number of kinetic species
+ADIAB_E = 0; % adiabatic electron model
+ADIAB_I = 1; % adiabatic ion model
+BETA = 1e-5; % electron plasma beta
+MHD_PD = 0; % MHD pressure drift
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 2; % real space x-gridpoints
+NY = 16; % real space y-gridpoints
+LX = 2*pi/0.3; % Size of the squared frequency domain in x direction
+LY = 2*pi/4.0; % Size of the squared frequency domain in y direction
+NZ = 16; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+% GEOMETRY= 'miller';
+GEOMETRY= 'z-pinch';
+EPS = 1.0; % inverse aspect ratio
+Q0 = 1.0; % safety factor
+SHEAR = 0.0; % magnetic shear
+KAPPA = 1.0; % elongation
+S_KAPPA = 0.0;
+DELTA = 0.0; % triangularity
+S_DELTA = 0.0;
+ZETA = 0.0; % squareness
+S_ZETA = 0.0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE= 0;
+%% TIME PARAMETERS
+TMAX = 10; % Maximal time unit
+DT = 1e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 3; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.1; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 10.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ % RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+ RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+
+if 1 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+
diff --git a/wk/KBM_ky_PJ_scan.m b/wk/KBM_ky_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..1fe7d44587a1e1b23b7d3e3b6cc7f25e8dca716e
--- /dev/null
+++ b/wk/KBM_ky_PJ_scan.m
@@ -0,0 +1,218 @@
+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 = 'gyacomo23_sp';
+% EXECNAME = 'gyacomo23_dp';
+% EXECNAME = 'gyacomo23_debug';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'lin_KBM'; % Name of the simulation
+RERUN = 1; % rerun if the data does not exist
+RUN = 1;
+K_Ne = 3; % ele Density '''
+K_Te = 4.5; % ele Temperature '''
+K_Ni = 3; % ion Density gradient drive
+K_Ti = 8; % ion Temperature '''
+P_a = [2 4 8 16];
+% P_a = 10;
+ky_a= 0.05:0.1:0.75;
+% ky_a = 0.05;
+% collision setting
+CO = 'DG';
+GKCO = 1; % gyrokinetic operator
+COLL_KCUT = 1.75;
+NU = 1e-2;
+% model
+NA = 2;
+BETA = 0.03; % electron plasma beta
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT0 = 5e-3;
+TMAX = 15;
+% arrays for the result
+g_ky = zeros(numel(ky_a),numel(P_a),2);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+j = 1;
+for P = P_a
+ J = P/2;
+i = 1;
+for ky = ky_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
+ DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = P/2; % Laguerre "
+ NX = 12; % real space x-gridpoints
+ NY = 2;
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/ky; % 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
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 0.5; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% 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 = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ ADIAB_E = (NA==1);
+ 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;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ Ntime = numel(data_.Ts0D);
+ catch
+ data_.outfilenames = [];
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
+ end
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>5
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % 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
+
+ [~,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 = 2: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(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(ky_a)>1 && numel(P_a)>1)
+%% Save metadata
+ pmin = num2str(min(P_a)); pmax = num2str(max(P_a));
+ kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
+filename = [num2str(NX),'x',num2str(NZ),'_ky_',kymin,'_',kymax,...
+ '_P_',pmin,'_',pmax,'_',CONAME,'_',num2str(NU),'_be_',num2str(BETA),'.mat'];
+metadata.name = filename;
+metadata.kymin = ky;
+metadata.title = ['$\nu_{',CONAME,'}=$',num2str(NU),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s2name = '$P$';
+metadata.s2 = P_a;
+metadata.s1name = '$ky$';
+metadata.s1 = ky_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/KBM_load_metadata_scan.m b/wk/KBM_load_metadata_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..5309527b6eff7ba71007b7e46860e9222313c21f
--- /dev/null
+++ b/wk/KBM_load_metadata_scan.m
@@ -0,0 +1,171 @@
+% Metadata path
+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';
+
+datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
+
+%% Chose if we filter gamma>0.05
+FILTERGAMMA = 0;
+
+%% Load data
+fname = ['../results/',datafname];
+d = load(fname);
+if FILTERGAMMA
+ d.data = d.data.*(d.data>0.025);
+ d.err = d.err.*(d.data>0.025);
+end
+if 0
+%% Pcolor of the peak
+figure;
+% [XX_,YY_] = meshgrid(d.s1,d.s2);
+[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
+pclr=imagesc_custom(XX_,YY_,d.data'.*(d.data>0)');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),d.data'.*(d.data>0)');
+title(d.title);
+xlabel(d.s1name); ylabel(d.s2name);
+set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
+set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
+colormap(jet)
+colormap(bluewhitered)
+clb=colorbar;
+clb.Label.String = '$\gamma c_s/R$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
+if 1
+%% Scan along first dimension
+figure
+colors_ = jet(numel(d.s2));
+for i = 1:numel(d.s2)
+ % plot(d.s1,d.data(:,i),'s-',...
+ plot(d.s1(d.data(:,i)>0),d.data((d.data(:,i)>0),i),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ 'color',colors_(i,:));
+ % errorbar(d.s1,d.data(:,i),d.err(:,i),'s-',...
+ % 'LineWidth',2.0,...
+ % 'DisplayName',[d.s2name,'=',num2str(d.s2(i))],...
+ % 'color',colors_(i,:));
+ hold on;
+end
+xlabel(d.s1name); ylabel(d.dname);title(d.title);
+xlim([d.s1(1) d.s1(end)]);
+colormap(colors_);
+clb = colorbar;
+% caxis([d.s2(1)-0.5,d.s2(end)+0.5]);
+clim([1 numel(d.s2)+1]);
+clb.Ticks=linspace(d.s2(1),d.s2(end),numel(d.s2));
+clb.Ticks =1.5:numel(d.s2)+1.5;
+clb.TickLabels=d.s2;
+clb.Label.String = d.s2name;
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
+if 0
+%% Scan along second dimension
+figure
+colors_ = jet(numel(d.s1));
+for i = 1:numel(d.s1)
+ plot(d.s2,d.data(i,:),'s-',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
+ 'color',colors_(i,:));
+% errorbar(d.s2,d.data(i,:),d.err(i,:),'s-',...
+% 'LineWidth',2.0,...
+% 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
+% 'color',colors_(i,:));
+ hold on;
+end
+xlabel(d.s2name); ylabel(d.dname);title(d.title);
+xlim([d.s2(1) d.s2(end)]);
+colormap(jet(numel(d.s1)));
+clb = colorbar;
+caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
+clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
+clb.YTick=d.s1;
+clb.Label.String = d.s1name;
+clb.TickLabelInterpreter = 'latex';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
+
+if 0
+%% Convergence analysis
+figure
+% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
+colors_ = jet(numel(d.s1));
+for i = 1:numel(d.s1)
+% target_ = d.data(i,end);
+ target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
+ if target_ > 0
+ eps_ = abs(target_ - d.data(i,1:end-1))/abs(target_);
+ semilogy(d.s2(1:end-1),eps_,'-s',...
+ 'LineWidth',2.0,...
+ 'DisplayName',[d.s1name,'=',num2str(d.s1(i))],...
+ 'color',colors_(i,:));
+ hold on;
+ end
+end
+xlabel(d.s2name); ylabel('$\epsilon_r$');title(d.title);
+xlim([d.s2(1) d.s2(end)]);
+colormap(colors_);
+clb = colorbar;
+caxis([d.s1(1)-0.5,d.s1(end)+0.5]);
+clb.Ticks=linspace(d.s1(1),d.s1(end),numel(d.s1));
+clb.YTick=d.s1;
+clb.Label.String = d.s1name;
+clb.TickLabelInterpreter = 'latex';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+grid on;
+end
+if 0
+%% Pcolor of the error
+figure; i_ = 0;
+% target_ = 2.72724991618068013377e-01; % Value for nuDGDK = 1.0, kT=6.96, (40,20), Nkx=8
+% target_ = 2.79666916212537142172e-01; % Value for nuDGDK = 0.05, kT=6.96, (40,20), Nkx=8
+% target_ = 2.71455e-01; % Value for nuDGDK = 0.001, kT=6.96, (40,20), Nkx=8
+% target_ = 1.39427e-01; % Value for nuDGDK = 0.001, kT=5.3, (30,16), Nkx=8
+target_ = 2.72510405826983714839e-01 % Value for nuDGDK = 0.001, kT=6.96, (50,25), Nkx=8
+% target_ = 2.73048910051283844069e-01; % Value for nuDGDK = 0.0, kT=6.96, (40,20), Nkx=8
+% target_ = 0.25*(d.data(end,end)+d.data(end-i_,end)+d.data(end,end-i_)+d.data(end-i_,end-i_));
+% eps_ = log(abs(target_ - d.data)/abs(target_));
+eps_ = max(-10,log(abs(target_ - d.data)/abs(target_)));
+sign_ = 1;%sign(d.data - target_);
+eps_ = d.data;
+for i = 1:numel(d.s1)
+ for j = 1:numel(d.s2)
+ % target_ = d.data(i,end);
+ % target_ = d.data(i,end);
+ % target_ = d.data(end,j);
+ eps_(i,j) = log(abs(target_ - d.data(i,j))/target);
+ if target_ > 0
+ % eps_(i,:) = max(-12,log(abs(target_ - d.data(i,1:end))/abs(target_)));
+ % eps_(i,:) = log(abs(target_ - d.data(i,1:end))/abs(target_));
+ % eps_(i,:) = min(100,100*abs(target_ - d.data(i,1:end))/abs(target_));
+ else
+ end
+ end
+end
+[XX_,YY_] = meshgrid(d.s1,d.s2);
+[XX_,YY_] = meshgrid(1:numel(d.s1),1:numel(d.s2));
+pclr=imagesc_custom(XX_,YY_,eps_'.*(d.data>0)'.*sign_');
+% pclr=contourf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_',5);
+% pclr=surf(1:numel(d.s1),1:numel(d.s2),eps_'.*(d.data>0)'.*sign_');
+title(d.title);
+xlabel(d.s1name); ylabel(d.s2name);
+set(gca,'XTick',1:numel(d.s1),'XTicklabel',d.s1)
+set(gca,'YTick',1:numel(d.s2),'YTicklabel',d.s2)
+% colormap(jet)
+colormap(bluewhitered)
+% caxis([-10, 0]);
+clb=colorbar;
+clb.Label.String = '$\log(\epsilon_r)$';
+clb.Label.Interpreter = 'latex';
+clb.Label.FontSize= 18;
+end
\ No newline at end of file
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index c36886926c94ac1ba6ad79bcc3df526042a34582..20777d70c099eddf7dd18ec7338862f53b312593 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -60,9 +60,9 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
% folder = '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_4.0/';
-folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
+% folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
% Miller CBC
-% folder = '/misc/gene_results/CBC/Miller_KT_6.96_128x64x24x16x8/';
+folder = '/misc/gene_results/Miller_KT_6.96_128x64x24x16x8/output/';
% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
% debug ? shearless
@@ -176,21 +176,21 @@ names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
figure;
subplot(311)
for i = 1:6
- plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
+ xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show'); title('GENE geometry');
subplot(312)
for i = 7:10
- plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+ xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show');
subplot(313)
for i = [11 12 14]
- plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ plot(gene_data.grids.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+ xlim([min(gene_data.grids.z),max(gene_data.grids.z)]); legend('show');
end
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 97ae221785047883d5dc3beeaedcbc7ec634cdf8..326ec67034e03bdd95d152cccb2fe11731fd8bd6 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -137,7 +137,7 @@ options.PLT_FTUBE = 0;
data.EPS = 0.4;
data.rho_o_R = 3e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
-save_figure(data,fig,'.png')
+% save_figure(data,fig,'.png')
end
if 0
@@ -231,7 +231,7 @@ end
if 0
%% Metric infos
-options.SHOW_FLUXSURF = 0;
+options.SHOW_FLUXSURF = 1;
options.SHOW_METRICS = 1;
[fig, geo_arrays] = plot_metric(data,options);
end
diff --git a/wk/benchmark_and_scan_scripts/Lin_1999_fig3_heat_diff_vs_time.txt b/wk/benchmark_and_scan_scripts/Lin_1999_fig3_heat_diff_vs_time.txt
new file mode 100644
index 0000000000000000000000000000000000000000..490ae85f1b52ee4e2ee534eb34858a770a68ff9b
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/Lin_1999_fig3_heat_diff_vs_time.txt
@@ -0,0 +1,85 @@
+22.939632545931715, 0.0028421861631937606
+46.876640419947535, 0.0015479763924537426
+73.80577427821521, 0.0016458578036862015
+91.75853018372703, 0.016904482243587093
+106.71916010498688, 0.043201954728034675
+117.69028871391077, 0.07639100361073647
+123.67454068241472, 0.11646800365423937
+129.65879265091866, 0.16897594292426155
+131.65354330708664, 0.218706950305245
+139.63254593175856, 0.32094589695625053
+143.62204724409452, 0.4217891271878308
+151.60104986876644, 0.5585584605791679
+155.5905511811024, 0.6994569394295327
+162.57217847769027, 0.8541784486883891
+169.5538057742782, 0.9453640463450356
+177.53280839895018, 0.7313046504546048
+182.51968503937013, 0.5172343788518148
+185.5118110236221, 0.4426596192050579
+196.482939632546, 0.365351430518699
+209.44881889763792, 0.2797631994895665
+225.4068241469817, 0.23700352373080436
+242.36220472440954, 0.21358448978408084
+256.3254593175854, 0.1597678397935065
+264.3044619422573, 0.12250402401357285
+285.249343832021, 0.09357462913820858
+308.1889763779529, 0.06050883832891052
+319.1601049868768, 0.06331114688012063
+330.13123359580055, 0.07301953277939699
+346.0892388451445, 0.06617145923057977
+357.0603674540684, 0.05239918214643069
+376.0104986876642, 0.05661170806687843
+390.97112860892406, 0.05804730209828746
+415.9055118110238, 0.05123185568654742
+450.81364829396335, 0.04721509258856482
+480.73490813648306, 0.05699235799944913
+503.67454068241494, 0.07917518597468143
+531.6010498687665, 0.08756398544104638
+553.5433070866145, 0.08902495613462691
+580.4724409448822, 0.11674714693812438
+600.4199475065618, 0.13477545279215797
+621.3648293963256, 0.14452009106596486
+636.3254593175855, 0.13214353040124127
+651.2860892388453, 0.16120343382491553
+661.2598425196852, 0.18748278012209796
+676.2204724409451, 0.21654268354577222
+687.191601049869, 0.27597482635112613
+690.1837270341209, 0.31189730427343
+694.1732283464569, 0.3464421919635735
+700.1574803149608, 0.3727070373109439
+702.1522309711288, 0.3768579341946898
+717.1128608923887, 0.3755310972868723
+723.0971128608926, 0.3548346166673918
+727.0866141732286, 0.3216999463464858
+736.0629921259845, 0.27062760110787254
+745.0393700787404, 0.24441713432229806
+758.0052493438322, 0.18507199721581769
+768.9763779527561, 0.1409129798001768
+777.952755905512, 0.10779643566653618
+785.931758530184, 0.08020112817389546
+792.9133858267719, 0.0636519192007079
+799.8950131233598, 0.0526275721059728
+820.8398950131236, 0.04856005568364741
+844.7769028871394, 0.04864706138252073
+855.7480314960633, 0.044543292585664584
+872.7034120734911, 0.047367352561592746
+897.6377952755906, 0.06403256913327837
+914.5931758530187, 0.07376270645727301
+929.5538057742783, 0.07934194689752161
+955.4855643044623, 0.08219863401052785
+969.4488188976381, 0.08639303374371021
+983.412073490814, 0.09611229535534571
+989.396325459318, 0.09751526224967744
+1008.3464566929135, 0.1031090036397384
+1027.2965879265096, 0.11422760690825251
+1037.2703412073495, 0.1253135830396891
+1047.2440944881894, 0.13087469729267276
+1061.2073490813652, 0.14888125172198774
+1068.1889763779532, 0.16686242948913155
+1076.167979002625, 0.1820848015545018
+1084.1469816272968, 0.19868838908948538
+1100.1049868766409, 0.22913313322022588
+1113.0708661417327, 0.24713606241208796
+1116.0629921259845, 0.24852815359406044
+1126.0367454068244, 0.2416583286205246
+1134.0157480314965, 0.2237315294151767
diff --git a/wk/benchmark_and_scan_scripts/old/fort_00.90 b/wk/benchmark_and_scan_scripts/old/fort_00.90
index b4c87cc3febeef70b10cb4afa9aeb4d00819323c..70572048a8317bae591be2e7b8413f7ad49218ce 100644
--- a/wk/benchmark_and_scan_scripts/old/fort_00.90
+++ b/wk/benchmark_and_scan_scripts/old/fort_00.90
@@ -1,23 +1,23 @@
&BASIC
nrun = 100000000
- dt = 0.02
- tmax = 30
+ dt = 0.005
+ tmax = 15
maxruntime = 356400
job2load = -1
/
&GRID
- pmax = 60
- jmax = 30
- Nx = 8
- Lx = 7.854
+ pmax = 4
+ jmax = 2
+ Nx = 9
+ Lx = 62.8319
Ny = 2
- Ly = 20.944
- Nz = 24
+ Ly = 25.1327
+ Nz = 16
SG = .false.
Nexc = 1
/
&GEOMETRY
- geom = 's-alpha'
+ geom = 'miller'
q0 = 1.4
shear = 0.8
eps = 0.18
@@ -34,31 +34,32 @@
dtsave_2d = -1
dtsave_3d = 2
dtsave_5d = 100
- write_doubleprecision = .false.
+ write_doubleprecision = .true.
write_gamma = .true.
write_hf = .true.
write_phi = .true.
write_Na00 = .true.
- write_Napj = .false.
- write_dens = .false.
+ write_Napj = .true.
+ write_dens = .true.
write_temp = .true.
/
&MODEL_PAR
- LINEARITY = 'linear'
- Na = 1
+LINEARITY = 'linear'
+RM_LD_T_EQ= .false.
+ Na = 2
mu_x = 0
mu_y = 0
N_HD = 4
- mu_z = 1
+ mu_z = 2
HYP_V = 'hypcoll'
mu_p = 0
mu_j = 0
- nu = 0.1
+ nu = 0.005
k_gB = 1
k_cB = 1
lambdaD = 0
- beta = 0.0001
- ADIAB_E = .true.
+ beta = 0.03
+ ADIAB_E = .false.
/
&CLOSURE_PAR
hierarchy_closure='truncation'
@@ -71,20 +72,27 @@
tau_ = 1
sigma_ = 1
q_ = 1
- K_N_ = 2.22
- K_T_ = 3.5
+ K_N_ = 3
+ K_T_ = 8
+/
+&SPECIES
+ name_ = electrons tau_ = 1
+ sigma_ = 0.051962
+ q_ = -1
+ K_N_ = 3
+ K_T_ = 4.5
/
&COLLISION_PAR
collision_model = 'DG'
GK_CO = .false.
INTERSPECIES = .true.
mat_file = '/home/ahoffman/gyacomo/wk/benchmark_and_scan_scripts/oiCa/null'
- collision_kcut = 1.75
+ collision_kcut = 1
/
&INITIAL_CON
INIT_OPT = 'phi'
- init_background = 0
- init_noiselvl = 1e-05
+ init_background = 1e-05
+ init_noiselvl = 0
iseed = 42
/
&TIME_INTEGRATION_PAR
diff --git a/wk/benchmark_and_scan_scripts/old/lin_KBM.m b/wk/benchmark_and_scan_scripts/old/lin_KBM.m
index 7a8c6b786c23f3f991675d41e4eee72bb20cec72..9d52a8a45d6e2d466ece1b28572afcbe16db0a33 100644
--- a/wk/benchmark_and_scan_scripts/old/lin_KBM.m
+++ b/wk/benchmark_and_scan_scripts/old/lin_KBM.m
@@ -1,188 +1,162 @@
%% 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 = 'test_circular_geom'; % Name of the simulation
-% SIMID = 'linear_CBC'; % Name of the simulation
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
SIMID = 'lin_KBM'; % Name of the simulation
-RUN = 1 ; % To run or just to load
-addpath(genpath('../matlab')) % ... add
+RUN = 1; % To run or just to load
default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-EXECNAME = 'helaz3';
-% EXECNAME = 'helaz3_dbg';
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% 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_Ne = 3; % ele Density '''
-K_Te = 4.5; % ele Temperature '''
-K_Ni = 3; % ion Density gradient drive
-K_Ti = 8; % ion Temperature '''
+EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.005; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 3; % ele Density '''
+K_Te = 4.5; % ele Temperature '''
+K_Ni = 3; % ion Density gradient drive
+K_Ti = 8; % 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.03; % electron plasma beta
-%% GRID PARAMETERS
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.03; % electron plasma beta
+%% Set up grid parameters
P = 4;
-J = P/2;
-PMAXE = P; % Hermite basis size of electrons
-JMAXE = J; % Laguerre "
-PMAXI = P; % " ions
-JMAXI = J; % "
-NX = 9; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
-LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.25; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
-NPOL = 1;
-SG = 0; % Staggered z grids option
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 9; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.25; % Size of the squared frequency domain in y direction
+NZ = 16; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
%% GEOMETRY
-% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
-GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
-NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
-EPS = 0.18; % inverse aspect ratio
-%% TIME PARMETERS
-TMAX = 15.1; % Maximal time unit
-DT = 5e-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/5; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
-JOB2LOAD= -1;
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 5e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
%% 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 = 0; % 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
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (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;
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% 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 = 10; %
-MU_P = 0.0; %
-MU_J = 0.0; %
-LAMBDAD = 0.0;
-NOISE0 = 0.0e-5; % Init noise amplitude
-BCKGD0 = 1.0; % Init background
-k_gB = 1.0;
-k_cB = 1.0;
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 2.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+
%%-------------------------------------------------------------------------
%% RUN
setup
% system(['rm fort*.90']);
% Run linear simulation
if RUN
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; time mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ',HELAZDIR,'bin/',EXECNAME,' 1 4 1 0; cd ../../../wk'])
- system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ',HELAZDIR,'bin/',EXECNAME,' 1 1 1 0; cd ../../../wk'])
-% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'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'])
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
end
-%% Load results
-%%
-filename = [SIMID,'/',PARAMS,'/'];
-LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
% 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
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
-%% Short analysis
if 0
-%% 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);
-[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.25];
-options.normalized = 1;
-% options.field = 'phi';
-fig = plot_ballooning(data,options);
+%% Plot heat flux evolution
+figure
+semilogy(data.Ts0D,data.HFLUX_X);
+xlabel('$tc_s/R$'); ylabel('$Q_x$');
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)
+if 1 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in 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 0
-%% RH TEST
-ikx = 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(1,ikx,:,it0:it1)),3))./squeeze(mean(real(x(1,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=0.00$',data.kx(ikx)))
-end
diff --git a/wk/benchmark_and_scan_scripts/plot_sparc_edge_profiles.m b/wk/benchmark_and_scan_scripts/plot_sparc_edge_profiles.m
new file mode 100644
index 0000000000000000000000000000000000000000..f086f7aa115956e8531066a7413d8a9d3bc49d00
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/plot_sparc_edge_profiles.m
@@ -0,0 +1,30 @@
+d_ = load("rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt");
+rho_ne = d_(:,1); ne = d_(:,2);
+d_ = load("rodriguez_fernandez_2020_fig3_sparc_T_prof.txt");
+rho_Te = d_(:,1); Te = d_(:,2);
+figure;
+plot(rho_ne,ne,'.k','DisplayName','$n_e$'); hold on;
+plot(rho_Te,Te,'or','DisplayName','$T_e$');
+%%
+[~,i1] = min(abs(rho_ne-0.95));
+[~,i2] = min(abs(rho_ne-0.98));
+[~,in975] = min(abs(rho_ne-0.975));
+
+fit_n = polyfit(rho_ne(i1:i2),ne(i1:i2),1);
+
+rho_ = 0.925:0.025:1.0;
+
+plot(rho_,fit_n(1).*rho_+fit_n(2),'-k','DisplayName',['$-\partial_x n/n_{975}=$',num2str(-fit_n(1)/ne(in975))]);
+
+%%
+[~,i1] = min(abs(rho_Te-0.95));
+[~,i2] = min(abs(rho_Te-0.98));
+[~,iT975] = min(abs(rho_Te-0.975));
+
+fit_T = polyfit(rho_Te(i1:i2),Te(i1:i2),1);
+
+rho_ = 0.925:0.025:1.0;
+
+plot(rho_,fit_T(1).*rho_+fit_T(2),'-r','DisplayName',['$-\partial_x T/T_{975}=$',num2str(-fit_T(1)/Te(iT975))]);
+xline(rho_Te(iT975));
+legend('show')
\ No newline at end of file
diff --git a/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_T_prof.txt b/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_T_prof.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fe16085b99efc5ab554e87935faaa302118ae522
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_T_prof.txt
@@ -0,0 +1,68 @@
+0.7368446165421814, 6.290726909030852
+0.7448753005083408, 6.144500926758546
+0.7468572791837185, 6.364353746490245
+0.75290231414362, 6.034904846671925
+0.758917986456479, 5.998495164339065
+0.768949000752418, 5.888972490870049
+0.7809840157090162, 5.779523224018643
+0.7870033583527556, 5.7064836395000995
+0.7950046796718726, 5.853296874713255
+0.795030371988035, 5.5968875594134815
+0.8030353636380321, 5.707070892440953
+0.8110880695894738, 5.3410654970545615
+0.8151070819034336, 5.231322603732725
+0.8250977225596889, 5.524728854306215
+0.8271237452056306, 5.305022847809731
+0.8351580995026704, 5.122166963351745
+0.837140078178048, 5.3420197830834475
+0.8411774421464093, 5.049127378833202
+0.8532197977647686, 4.866418307610431
+0.8552164577636674, 4.939751518599397
+0.8592317997467473, 4.866638527463252
+0.8632398010680664, 4.866785340698463
+0.8692554733809255, 4.830375658365604
+0.8913251729643428, 4.574773815859498
+0.9073681892422603, 4.465471362243306
+0.9113761905635795, 4.465618175478522
+0.9153731808922576, 4.575654695270778
+0.915387862215779, 4.429135086528049
+0.9294342184947975, 4.246499421922888
+0.9294342184947975, 4.246499421922888
+0.9354425501458956, 4.283349543961386
+0.945477234772715, 4.137196968306689
+0.9494925767557947, 4.064083977170544
+0.953533611055037, 3.7345616707346174
+0.957508579398433, 4.0643776036409704
+0.957508579398433, 4.0643776036409704
+0.957552623368997, 3.6248187774127842
+0.9595309317134941, 3.881301499330167
+0.9595676350222972, 3.5150024774733453
+0.9596043383311006, 3.1487034556165234
+0.9596116789928613, 3.075443651245159
+0.9615642950211962, 3.5883356884623154
+0.961652382962324, 2.7092180360059466
+0.9635866473362575, 3.405259584151512
+0.9636784056082659, 2.4895120295094593
+0.9656163403130793, 3.1489236754693444
+0.9656493732910023, 2.8192545557982065
+0.9656970875924467, 2.34306582738434
+0.9677194399075078, 2.1599897230735365
+0.969686737259364, 2.526362151547964
+0.9697528032152098, 1.8670239122056849
+0.97165036428034, 2.9293644822080687
+0.9737094199042047, 2.379989356040447
+0.9737975078453325, 1.5008717035840782
+0.9757354425501461, 2.1602833495439633
+0.9758125194986329, 1.3910554036446392
+0.977757794865207, 1.9772072452331564
+0.9777724761887283, 1.8306876364904312
+0.9818024994953296, 1.6110550366115497
+0.9818355324732526, 1.2813859169404118
+0.9818648951202953, 0.9883466994549543
+0.9838578847883137, 1.098309812629605
+0.9879062597493165, 0.6955277018223178
+0.9879062597493165, 0.6955277018223178
+0.9898882384246943, 0.9153805215540203
+0.989895579086455, 0.8421207171826559
+0.9939182617312952, 0.6957479216751388
+0.9999449450367952, 0.549448532785231
diff --git a/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt b/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e148577b821f81d4c81cbd8fde8597b7611840fa
--- /dev/null
+++ b/wk/benchmark_and_scan_scripts/rodriguez_fernandez_2020_fig3_sparc_ne_prof.txt
@@ -0,0 +1,43 @@
+0.8907345182249297, 2.9380728997180983
+0.9087345840386986, 2.9383361547939852
+0.9147309496561232, 2.9237986244922283
+0.9227236662656902, 2.894665062760742
+0.9307236955162542, 2.8947820650166918
+0.9367200611336789, 2.8802445347149357
+0.9407164194384622, 2.8656777538491927
+0.9447091214227475, 2.83648569098973
+0.9486945107660357, 2.778043064142831
+0.9526616185068318, 2.64647402732734
+0.9526799001093241, 2.7196004372959317
+0.9526799001093241, 2.7196004372959317
+0.95664700785012, 2.588031400480441
+0.9585994829962815, 2.397931985126089
+0.9586323898807673, 2.5295595230695547
+0.9605885213474273, 2.3540853897089207
+0.9625592780960809, 2.2371123843231615
+0.962570247057576, 2.280988230304315
+0.9625848723395699, 2.33948935827919
+0.9645263785242358, 2.1055140969436827
+0.9645300348447342, 2.1201393789374023
+0.9645519727677249, 2.207891070899713
+0.9645519727677249, 2.207891070899713
+0.9665117605548833, 2.047042219532795
+0.9684825173035368, 1.9300692141470357
+0.9685007989060289, 2.0031956241156283
+0.9704569303726888, 1.8277214907549926
+0.9723509043908755, 1.4036175635011467
+0.9723655296728692, 1.4621186914760216
+0.9723801549548629, 1.5206198194508946
+0.9723911239163582, 1.5644956654320499
+0.9724057491983519, 1.6229967934069247
+0.9724203744803456, 1.6814979213817978
+0.9724386560828379, 1.7546243313503904
+0.9763106994906746, 1.2427979626982193
+0.9763289810931667, 1.3159243726668102
+0.9802741509109724, 1.0966036438890097
+0.9802851198724677, 1.140479489870165
+0.9822558766211212, 1.023506484484404
+0.9842485712927652, 0.9942851710609553
+0.9842485712927652, 0.9942851710609553
+0.9882266479950569, 0.9065919802266187
+0.9882303043155553, 0.9212172622203383
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index e73ed755311964f45701a830945a5f9ece53b3f6..0c58560405fa2baf611637c4d187dd78a91cf186 100644
--- a/wk/fast_analysis.m
+++ b/wk/fast_analysis.m
@@ -5,8 +5,8 @@ addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
addpath(genpath([gyacomodir,'matlab/load'])) % ... add
default_plots_options
% Partition of the computer where the data have to be searched
-PARTITION = '/misc/gyacomo23_outputs/';
-% PARTITION = '/home/ahoffman/gyacomo/';
+% PARTITION = '/misc/gyacomo23_outputs/';
+PARTITION = '/home/ahoffman/gyacomo/';
%% Tests
% resdir = 'test_stepon_AA/CBC_stepon_AA';
@@ -56,11 +56,21 @@ PARTITION = '/misc/gyacomo23_outputs/';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.2';
-% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/9x5x192x96x24/nu_0.02';
-resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x192x96x24/nu_0.2';
+
+% resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/continue_SG_nu_0.02';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/5x3x192x96x24/nu_0.02';
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/9x5x192x96x24/nu_0.02';
% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/17x9x192x96x24/nu_0.02';
+
+
+% resdir = 'paper_2_GYAC23/collision_study/nuSGGK_scan_kT_5.3/17x9x192x96x24/nu_0.005';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/5x3x192x96x24/nu_0.005';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/9x5x192x96x24/nu_0.005';
+% resdir = 'paper_2_GYAC23/collision_study/nuLDGK_scan_kT_5.3/17x9x128x64x24/nu_0.005';
%% kT eff study
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx120';
% resdir = 'paper_2_GYAC23/kT_eff_study/1x3x128x64x24_kT_3.0/Lx240';
@@ -71,9 +81,9 @@ resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.
% resdir = 'paper_2_GYAC23/kT_eff_study/7x3x128x64x24_kT_3.6/dmax';
%% dev
-% PARTITION='';
+PARTITION='';
% resdir = '/home/ahoffman/gyacomo/testcases/zpinch_3D';
-% resdir = '/home/ahoffman/gyacomo';
+resdir = '/home/ahoffman/gyacomo/results/dev/3D_kine_zpinch_test';
%% CBC benchmark
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/3x2x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/5x3x128x64x24/kT_7.0';
@@ -81,7 +91,7 @@ resdir = 'paper_2_GYAC23/collision_study/nuDGGK_scan_kT_5.3/17x9x192x96x24/nu_0.
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/17x9x128x64x24/kT_7.0';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/31x16x128x64x24/kT_7.0';
%%
-J0 = 00; J1 = 10;
+J0 = 00; J1 = 20;
% Load basic info (grids and time traces)
DATADIR = [PARTITION,resdir,'/'];
@@ -111,16 +121,16 @@ if 0
% data.Ni00 = reshape(data.Na00(1,:,:,:,:),data.grids.Nky,data.grids.Nkx,data.grids.Nz,numel(data.Ts3D));
options.INTERP = 1;
options.POLARPLOT = 0;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = '\phi^{NZ}';
% options.NAME = '\omega_z';
-options.NAME = 'N_i^{00}';
+% options.NAME = 'N_i^{00}';
% options.NAME = 's_{Ey}';
% options.NAME = 'n_i^{NZ}';
% options.NAME = 'Q_x';
% options.NAME = 'n_i';
% options.NAME = 'n_i-n_e';
-options.PLAN = 'xy';
+options.PLAN = 'xz';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
@@ -144,8 +154,8 @@ options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 0;
options.NORMALIZE = 0;
-options.NAME = 'N_i^{00}';
-% options.NAME = '\phi';
+% options.NAME = 'N_i^{00}';
+options.NAME = '\phi';
options.PLAN = 'xy';
options.COMP = 'avg';
options.TIME = [100 200 300];
@@ -166,8 +176,8 @@ options.ST = 1;
options.NORMALIZED = 0;
options.LOGSCALE = 1;
options.FILTER = 0; %filter the 50% time-average of the spectrum from
-options.TAVG_2D = 1; %Show a 2D plot of the modes, 50% time averaged
-options.TAVG_2D_CTR= 1; %make it contour plot
+options.TAVG_2D = 0; %Show a 2D plot of the modes, 50% time averaged
+options.TAVG_2D_CTR= 0; %make it contour plot
fig = show_moments_spectrum(data,options);
end
diff --git a/wk/lin_AUG.m b/wk/lin_AUG.m
new file mode 100644
index 0000000000000000000000000000000000000000..ae21fb232f1cdfe21f58f263c1e5b65a76418c0e
--- /dev/null
+++ b/wk/lin_AUG.m
@@ -0,0 +1,177 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'lin_KBM'; % Name of the simulation
+RUN = 1; % To run or just to load
+default_plots_options
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
+% EXECNAME = 'gyacomo23_debug'; % single precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.1; % Collision frequency
+TAU = 1.38; % e/i temperature ratio
+K_Ne = 34; % ele Density '''
+K_Te = 56; % ele Temperature '''
+K_Ni = 34; % ion Density gradient drive
+K_Ti = 21; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 1.52e-4; % electron plasma beta
+MHD_PD = 0; % MHD pressure drift
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 4; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.3; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.2; % Size of the squared frequency domain in y direction
+NZ = 16; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+EPS = 0.3; % inverse aspect ratio
+% EPS = 0.18; % inverse aspect ratio
+Q0 = 5.7; % safety factor
+% Q0 = 1.4; % safety factor
+SHEAR = 4.6; % magnetic shear
+% SHEAR = 0.8; % magnetic shear
+KAPPA = 1.8; % elongation
+% KAPPA = 1.0; % elongation
+S_KAPPA = 0.0;
+DELTA = 0.4; % triangularity
+% DELTA = 0.0; % triangularity
+S_DELTA = 0.0;
+ZETA = 0.0; % squareness
+S_ZETA = 0.0;
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+PB_PHASE= 0;
+%% TIME PARAMETERS
+TMAX = 5; % Maximal time unit
+DT = 2e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 3; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'phi'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.1; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 10.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ % RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+ RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+
+if 0 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+if 1
+%% Ballooning plot
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+if data.inputs.BETA > 0
+[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
+end
+options.time_2_plot = [120];
+options.kymodes = [0.5];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
diff --git a/wk/lin_ITG_salpha.m b/wk/lin_ITG_salpha.m
index c603b379acbdeb9be97143e08806feabba42b9cc..8616c009349c73bd3b1fdb1b4a430321a350b411 100644
--- a/wk/lin_ITG_salpha.m
+++ b/wk/lin_ITG_salpha.m
@@ -20,19 +20,19 @@ EXECNAME = 'gyacomo23_sp'; % single precision
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-NU = 0.001; % Collision frequency
+NU = 0.005; % Collision frequency
TAU = 1.0; % e/i temperature ratio
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 2.22; % ion Density gradient drive
-K_Ti = 6.96; % ion Temperature
+K_Ti = 5.3; % ion Temperature
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
BETA = 0.0; % electron plasma beta
%% Set up grid parameters
-P = 2;
-J = 1;%P/2;
+P = 4;
+J = 2;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 8; % real space x-gridpoints
@@ -58,7 +58,7 @@ NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
TMAX = 50; % Maximal time unit
-DT = 10e-3; % Time step
+DT = 1e-2; % Time step
DTSAVE0D = 1; % Sampling per time unit for 0D arrays
DTSAVE2D = -1; % Sampling per time unit for 2D arrays
DTSAVE3D = 2; % Sampling per time unit for 3D arrays
@@ -68,7 +68,7 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-GKCO = 0; % Gyrokinetic operator
+GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
@@ -107,7 +107,7 @@ NOISE0 = 0.0e-5; % Initial noise amplitude
BCKGD0 = 1.0e-5; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
-COLL_KCUT = 1000; % Cutoff for collision operator
+COLL_KCUT = 1; % Cutoff for collision operator
%%-------------------------------------------------------------------------
%% RUN
diff --git a/wk/lin_KBM.m b/wk/lin_KBM.m
new file mode 100644
index 0000000000000000000000000000000000000000..40b3db95a90a41493ddc721ff55cc747b698e9df
--- /dev/null
+++ b/wk/lin_KBM.m
@@ -0,0 +1,168 @@
+%% QUICK RUN SCRIPT
+% This script creates a directory in /results and runs a simulation directly
+% from the Matlab framework. It is meant to run only small problems in linear
+% for benchmarking and debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab module
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot module
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute module
+addpath(genpath([gyacomodir,'matlab/load'])) % Add load module
+
+%% Set simulation parameters
+SIMID = 'lin_KBM'; % Name of the simulation
+RUN = 1; % To run or just to load
+default_plots_options
+EXECNAME = 'gyacomo23_sp'; % single precision
+% EXECNAME = 'gyacomo23_dp'; % double precision
+
+%% Set up physical parameters
+CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
+NU = 0.00; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 3; % ele Density '''
+K_Te = 4.5; % ele Temperature '''
+K_Ni = 3; % ion Density gradient drive
+K_Ti = 8; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = (NA==1); % adiabatic electron model
+BETA = 0.03; % electron plasma beta
+%% Set up grid parameters
+P = 2;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 12; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
+LY = 2*pi/0.25; % Size of the squared frequency domain in y direction
+NZ = 24; % number of perpendicular planes (parallel grid)
+SG = 0; % Staggered z grids option
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+
+%% GEOMETRY
+% GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+KAPPA = 1.0; % elongation
+DELTA = 0.0; % triangularity
+ZETA = 0.0; % squareness
+PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
+SHIFT_Y = 0.0; % Shift in the periodic BC in z
+NPOL = 1; % Number of poloidal turns
+
+%% TIME PARAMETERS
+TMAX = 15; % Maximal time unit
+DT = 5e-3; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+JOB2LOAD = -1; % Start a new simulation serie
+
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 0; % Gyrokinetic operator
+ABCO = 1; % INTERSPECIES collisions
+INIT_ZF = 0; % Initialize zero-field quantities
+HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+DMAX = -1;
+NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+NMAX = 0;
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT = 'mom00'; % Start simulation with a noisy mom00/phi/allmom
+NUMERICAL_SCHEME = 'RK4'; % Numerical integration scheme (RK2,SSPx_RK2,RK3,SSP_RK3,SSPx_RK3,IMEX_SSP2,ARK2,RK4,DOPRI5)
+
+%% OUTPUTS
+W_DOUBLE = 1; % Output flag for double moments
+W_GAMMA = 1; % Output flag for gamma (Gyrokinetic Energy)
+W_HF = 1; % Output flag for high-frequency potential energy
+W_PHI = 1; % Output flag for potential
+W_NA00 = 1; % Output flag for nalpha00 (density of species alpha)
+W_DENS = 1; % Output flag for total density
+W_TEMP = 1; % Output flag for temperature
+W_NAPJ = 1; % Output flag for nalphaparallel (parallel momentum of species alpha)
+W_SAPJ = 0; % Output flag for saparallel (parallel current of species alpha)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% UNUSED PARAMETERS
+% These parameters are usually not to play with in linear runs
+MU = 0.0; % Hyperdiffusivity coefficient
+MU_X = MU; % Hyperdiffusivity coefficient in x direction
+MU_Y = MU; % Hyperdiffusivity coefficient in y direction
+N_HD = 4; % Degree of spatial-hyperdiffusivity
+MU_Z = 1.0; % Hyperdiffusivity coefficient in z direction
+HYP_V = 'hypcoll'; % Kinetic-hyperdiffusivity model
+MU_P = 0.0; % Hyperdiffusivity coefficient for Hermite
+MU_J = 0.0; % Hyperdiffusivity coefficient for Laguerre
+LAMBDAD = 0.0; % Lambda Debye
+NOISE0 = 0.0e-5; % Initial noise amplitude
+BCKGD0 = 1.0e-5; % Initial background
+k_gB = 1.0; % Magnetic gradient strength
+k_cB = 1.0; % Magnetic curvature strength
+COLL_KCUT = 1; % Cutoff for collision operator
+
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
+% RUN =['time mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
+ RUN =['time mpirun -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+% RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 6 1 0;'];
+% RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ MVOUT='cd ../../../wk;';
+ system([MVIN,RUN,MVOUT]);
+end
+
+%% Analysis
+% load
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+
+if 0 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KX_TW = [0.2 1]*data.Ts3D(end);
+options.KY_TW = [0.2 1]*data.Ts3D(end);
+options.NMA = 1; % Set NMA option to 1
+options.NMODES = 999; % Set how much modes we study
+options.iz = 'avg'; % Compressing z
+options.ik = 1; %
+options.fftz.flag = 0; % Set fftz.flag option to 0
+fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+end
+
+if 1
+%% Ballooning plot
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+if data.inputs.BETA > 0
+[data.PSI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'psi');
+end
+options.time_2_plot = [120];
+options.kymodes = [0.25];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
diff --git a/wk/CBC_P_J_scan.m b/wk/paper_2_scripts_and_results/CBC_P_J_scan.m
similarity index 100%
rename from wk/CBC_P_J_scan.m
rename to wk/paper_2_scripts_and_results/CBC_P_J_scan.m
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/paper_2_scripts_and_results/CBC_kT_PJ_scan.m
similarity index 100%
rename from wk/CBC_kT_PJ_scan.m
rename to wk/paper_2_scripts_and_results/CBC_kT_PJ_scan.m
diff --git a/wk/CBC_kT_nu_scan.m b/wk/paper_2_scripts_and_results/CBC_kT_nu_scan.m
similarity index 100%
rename from wk/CBC_kT_nu_scan.m
rename to wk/paper_2_scripts_and_results/CBC_kT_nu_scan.m
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/paper_2_scripts_and_results/CBC_ky_PJ_scan.m
similarity index 99%
rename from wk/CBC_ky_PJ_scan.m
rename to wk/paper_2_scripts_and_results/CBC_ky_PJ_scan.m
index 0d3a80980a69016ea6421f0341d2c6360f150dbe..766104d7af3e366803c66a496f1e42afea554812 100644
--- a/wk/CBC_ky_PJ_scan.m
+++ b/wk/paper_2_scripts_and_results/CBC_ky_PJ_scan.m
@@ -11,14 +11,14 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
-RUN = 1;
+RUN = 0;
K_T = 5.3;
-P_a = [2 4 8 10];
+P_a = [2 4 8 16];
% P_a = 10;
ky_a= [0.05:0.05:1.0];
% ky_a = 0.05;
% collision setting
-CO = 'LD';
+CO = 'SG';
GKCO = 1; % gyrokinetic operator
COLL_KCUT = 1.75;
NU = 1e-2;
diff --git a/wk/CBC_ky_PJ_scan_colless.m b/wk/paper_2_scripts_and_results/CBC_ky_PJ_scan_colless.m
similarity index 96%
rename from wk/CBC_ky_PJ_scan_colless.m
rename to wk/paper_2_scripts_and_results/CBC_ky_PJ_scan_colless.m
index 37779f3fc8797a637df89fc069cd2f888ecdc5af..bb0f6c7f73b7be0fbe1c7858d69b887843692c67 100644
--- a/wk/CBC_ky_PJ_scan_colless.m
+++ b/wk/paper_2_scripts_and_results/CBC_ky_PJ_scan_colless.m
@@ -11,12 +11,12 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%
SIMID = 'p2_linear_new'; % Name of the simulation
RERUN = 0; % rerun if the data does not exist
-RUN = 1;
-K_T = 6.96;
-% P_a = [2 4 8 16 32 48];
-P_a = 64;
-% ky_a= [0.05:0.05:1.0];
-ky_a= [0.25:0.05:0.65];
+RUN = 0;
+K_T = 5.3;
+P_a = [2 4 8 16 32 64];
+% P_a = 32;
+ky_a= [0.05:0.05:1.0];
+% ky_a= [0.25:0.05:0.65];
% collision setting
CO = 'DG';
GKCO = 0; % gyrokinetic operator
@@ -32,10 +32,10 @@ K_N = 2.22; % Density '''
GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
-DT0 = 2e-3;
+DT0 = 1e-3;
TMAX = 30;
% arrays for the result
-g_ky = zeros(numel(ky_a),numel(P_a),NY/2+1);
+g_ky = zeros(numel(ky_a),numel(P_a),1);
g_avg= g_ky*0;
g_std= g_ky*0;
% Naming of the collision operator
@@ -46,6 +46,7 @@ else
end
j = 1;
for P = P_a
+J = P/2;
i = 1;
for ky = ky_a
%% PHYSICAL PARAMETERS
@@ -139,8 +140,8 @@ for ky = ky_a
end
if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
% system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 0; cd ../../../wk'])
end
data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
[data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
diff --git a/wk/paper_2_scripts_and_results/CBC_ky_nu_scan_colless.m b/wk/paper_2_scripts_and_results/CBC_ky_nu_scan_colless.m
new file mode 100644
index 0000000000000000000000000000000000000000..126ec7bd60511e62c139c03a101ab181c30aefe9
--- /dev/null
+++ b/wk/paper_2_scripts_and_results/CBC_ky_nu_scan_colless.m
@@ -0,0 +1,229 @@
+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 = 'gyacomo23_sp';
+% EXECNAME = 'gyacomo23_dp';
+% EXECNAME = 'gyacomo23_debug';
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%
+SIMID = 'p2_linear_new'; % Name of the simulation
+RERUN = 0; % rerun if the data does not exist
+RUN = 1;
+K_T = 5.3;
+P = 16;
+nu_a = [0.005 0.01 0.02 0.05];
+% nu_a = 0.05;
+ky_a = [0.05:0.05:1.0];
+% ky_a = 0.35:0.05:0.6;
+% collision setting
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
+COLL_KCUT = 1.00;
+% model
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 1e-4; % electron plasma beta
+% background gradients setting
+K_N = 2.22; % Density '''
+% Geometry
+% GEOMETRY= 'miller';
+GEOMETRY= 's-alpha';
+SHEAR = 0.8; % magnetic shear
+% time and numerical grid
+DT = 1e-3;
+TMAX = 40;
+% arrays for the result
+g_ky = zeros(numel(ky_a),numel(nu_a),1);
+g_avg= g_ky*0;
+g_std= g_ky*0;
+% Naming of the collision operator
+if GKCO
+ CONAME = [CO,'GK'];
+else
+ CONAME = [CO,'DK'];
+end
+j = 1;
+for NU = nu_a
+J = P/2;
+i = 1;
+for ky = ky_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)
+ K_Te = K_T; % ele Temperature '''
+ K_Ti = K_T; % ion Temperature '''
+ K_Ne = K_N; % ele Density '''
+ K_Ni = K_N; % ion Density gradient drive
+ %% GRID PARAMETERS
+ % DT = DT0/sqrt(J);
+ PMAX = P; % Hermite basis size
+ JMAX = P/2; % Laguerre "
+ NX = 8; % real space x-gridpoints
+ NY = 2;
+ LX = 2*pi/0.8; % Size of the squared frequency domain
+ LY = 2*pi/ky; % 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
+ DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+ DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+ DTSAVE3D = 2; % Sampling per time unit for 3D arrays
+ DTSAVE5D = 100; % Sampling per time unit for 5D arrays
+ JOB2LOAD = -1; % Start a new simulation serie
+ %% 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 = 0;
+ W_GAMMA = 1; W_HF = 1;
+ W_PHI = 1; W_NA00 = 1;
+ W_DENS = 0; W_TEMP = 1;
+ W_NAPJ = 0; W_SAPJ = 0;
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % unused
+ NA = 1;
+ ADIAB_E = (NA==1);
+ 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;
+ HYP_V = 'none';
+ HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+ DMAX = -1;
+ NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
+ NMAX = 0;
+ MU_Z = 1.0; %
+ MU_P = 0.0; %
+ MU_J = 0.0; %
+ LAMBDAD = 0.0;
+ NOISE0 = 1.0e-5; % Init noise amplitude
+ BCKGD0 = 0.0; % Init background
+ k_gB = 1.0;
+ k_cB = 1.0;
+ %% RUN
+ setup
+ % naming
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+ EXIST= 0;
+ % check if data exist to run if no data
+ data_ = {};
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ Ntime = numel(data_.Ts0D);
+ EXIST = 1;
+ catch
+ data_.outfilenames = [];
+ EXIST = 0;
+ end
+ if RUN && (RERUN || isempty(data_.outfilenames) || Ntime < 10)
+ % system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 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,' 3 2 1 0; cd ../../../wk'])
+ end
+ try
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+ if numel(data_.Ts3D)>10
+ if numel(data_.Ts3D)>10
+ % Load results after trying to run
+ filename = [SIMID,'/',PARAMS,'/'];
+ LOCALDIR = [gyacomodir,'results/',filename,'/'];
+
+ data_ = compile_results_low_mem(data_,LOCALDIR,00,00);
+ [data_.PHI, data_.Ts3D] = compile_results_3D(LOCALDIR,00,00,'phi');
+
+ % 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
+
+ [~,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 = 2: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(it1:it2));
+ tw = double(data_.Ts3D(it1:it2));
+ % gr = polyfit(tw,to_measure,1);
+ gr = fit(tw,to_measure,'poly1');
+ err= confint(gr);
+ g_ky(i,j,ik) = gr.p1;
+ g_std(i,j,ik) = abs(err(2,1)-err(1,1))/2;
+ end
+ [gmax, ikmax] = max(g_ky(i,j,:));
+
+ msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,data_.grids.ky(ikmax)); disp(msg);
+ end
+ end
+ catch
+ g_ky(i,j,:) = gr.p1;
+ g_std(i,j,:) = abs(err(2,1)-err(1,1))/2;
+ end
+ i = i + 1;
+end
+j = j + 1;
+end
+
+if 0
+%% Check time evolution
+figure;
+to_measure = log(field_t);
+plot(data_.Ts3D,to_measure); hold on
+plot(data_.Ts3D(it1:it2),to_measure(it1:it2),'--');
+end
+
+%% take max growth rate among z coordinate
+y_ = g_ky(:,:,2);
+e_ = g_std(:,:,2);
+
+%%
+if(numel(ky_a)>1 && numel(nu_a)>1)
+%% Save metadata
+ numin = num2str(min(nu_a)); numax = num2str(max(nu_a));
+ kymin = num2str(min(ky_a)); kymax= num2str(max(ky_a));
+filename = [num2str(NX),'x',num2str(NZ),'_P_',num2str(P),'_ky_',kymin,'_',kymax,...
+ '_',CONAME,'_nu_',numin,'_',numax,'_kT_',num2str(K_T),'.mat'];
+metadata.name = filename;
+metadata.kymin = ky;
+metadata.title = ['$P=$',num2str(P),'$\kappa_T=$',num2str(K_T),', $\kappa_N=$',num2str(K_N)];
+metadata.par = [num2str(NX),'x1x',num2str(NZ)];
+metadata.nscan = 2;
+metadata.s2name = ['$\nu_{',CONAME,'}=$'];
+metadata.s2 = nu_a;
+metadata.s1name = '$k_y$';
+metadata.s1 = ky_a;
+metadata.dname = '$\gamma c_s/R$';
+metadata.data = y_;
+metadata.err = e_;
+save([SIMDIR,filename],'-struct','metadata');
+disp(['saved in ',SIMDIR,filename]);
+clear metadata tosave
+end
diff --git a/wk/CBC_nu_PJ_scan.m b/wk/paper_2_scripts_and_results/CBC_nu_PJ_scan.m
similarity index 100%
rename from wk/CBC_nu_PJ_scan.m
rename to wk/paper_2_scripts_and_results/CBC_nu_PJ_scan.m
diff --git a/wk/Dimits_fig3.m b/wk/paper_2_scripts_and_results/Dimits_fig3.m
similarity index 100%
rename from wk/Dimits_fig3.m
rename to wk/paper_2_scripts_and_results/Dimits_fig3.m
diff --git a/wk/fast_heat_flux_analysis.m b/wk/paper_2_scripts_and_results/fast_heat_flux_analysis.m
similarity index 100%
rename from wk/fast_heat_flux_analysis.m
rename to wk/paper_2_scripts_and_results/fast_heat_flux_analysis.m
diff --git a/wk/paper_2_scripts_and_results/fort_00.90 b/wk/paper_2_scripts_and_results/fort_00.90
new file mode 100644
index 0000000000000000000000000000000000000000..c5c6dcd5e73b155d96dc82f3224793a73586753e
--- /dev/null
+++ b/wk/paper_2_scripts_and_results/fort_00.90
@@ -0,0 +1,93 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.0017678
+ tmax = 50
+ maxruntime = 356400
+ job2load = -1
+/
+&GRID
+ pmax = 16
+ jmax = 8
+ Nx = 8
+ Lx = 7.854
+ Ny = 2
+ Ly = 6.2832
+ Nz = 24
+ SG = .false.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 's-alpha'
+ q0 = 1.4
+ shear = 0.8
+ eps = 0.18
+ kappa = 1
+ delta = 0
+ zeta = 0
+ parallel_bc = 'dirichlet'
+ shift_y = 0
+ Npol = 1
+/
+&OUTPUT_PAR
+ dtsave_0d = 1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 2
+ dtsave_5d = 100
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .false.
+ write_dens = .false.
+ write_temp = .true.
+/
+&MODEL_PAR
+LINEARITY = 'linear'
+RM_LD_T_EQ= .false.
+ Na = 1
+ mu_x = 0
+ mu_y = 0
+ N_HD = 4
+ mu_z = 1
+ HYP_V = 'none'
+ mu_p = 0
+ mu_j = 0
+ nu = 0.01
+ k_gB = 1
+ k_cB = 1
+ lambdaD = 0
+ beta = 0.0001
+ ADIAB_E = .true.
+/
+&CLOSURE_PAR
+ hierarchy_closure='truncation'
+ dmax =-1
+ nonlinear_closure='truncation'
+ nmax =0
+/
+&SPECIES
+ name_ = ions
+ tau_ = 1
+ sigma_ = 1
+ q_ = 1
+ K_N_ = 2.22
+ K_T_ = 5.3
+/
+&COLLISION_PAR
+ collision_model = 'SG'
+ GK_CO = .true.
+ INTERSPECIES = .true.
+ mat_file = '/home/ahoffman/gyacomo/wk/paper_2_scripts_and_resuliCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+ collision_kcut = 1.75
+/
+&INITIAL_CON
+ INIT_OPT = 'phi'
+ init_background = 0
+ init_noiselvl = 1e-05
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/paper_2_scripts_and_results/header_3D_results.m
similarity index 100%
rename from wk/header_3D_results.m
rename to wk/paper_2_scripts_and_results/header_3D_results.m
diff --git a/wk/heat_flux_convergence_analysis.m b/wk/paper_2_scripts_and_results/heat_flux_convergence_analysis.m
similarity index 100%
rename from wk/heat_flux_convergence_analysis.m
rename to wk/paper_2_scripts_and_results/heat_flux_convergence_analysis.m
diff --git a/wk/kTcrit_convergence_analysis.m b/wk/paper_2_scripts_and_results/kTcrit_convergence_analysis.m
similarity index 100%
rename from wk/kTcrit_convergence_analysis.m
rename to wk/paper_2_scripts_and_results/kTcrit_convergence_analysis.m
diff --git a/wk/multiple_low_res_heat_flux_p2_analysis.m b/wk/paper_2_scripts_and_results/multiple_low_res_heat_flux_p2_analysis.m
similarity index 100%
rename from wk/multiple_low_res_heat_flux_p2_analysis.m
rename to wk/paper_2_scripts_and_results/multiple_low_res_heat_flux_p2_analysis.m
diff --git a/wk/multiple_sim_analysis.m b/wk/paper_2_scripts_and_results/multiple_sim_analysis.m
similarity index 100%
rename from wk/multiple_sim_analysis.m
rename to wk/paper_2_scripts_and_results/multiple_sim_analysis.m
diff --git a/wk/new_CBC_convergence_GENE_GYAC_comparison.m b/wk/paper_2_scripts_and_results/new_CBC_convergence_GENE_GYAC_comparison.m
similarity index 100%
rename from wk/new_CBC_convergence_GENE_GYAC_comparison.m
rename to wk/paper_2_scripts_and_results/new_CBC_convergence_GENE_GYAC_comparison.m
diff --git a/wk/nonlin_kT_scan_analysis.m b/wk/paper_2_scripts_and_results/nonlin_kT_scan_analysis.m
similarity index 100%
rename from wk/nonlin_kT_scan_analysis.m
rename to wk/paper_2_scripts_and_results/nonlin_kT_scan_analysis.m
diff --git a/wk/nonlin_nu_scan_analysis.m b/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
similarity index 94%
rename from wk/nonlin_nu_scan_analysis.m
rename to wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
index 8b64ff2efa9661837b8a5e38ead4633ff37050d2..60627f79cbdf77bd3abe8fd568748618087381d6 100644
--- a/wk/nonlin_nu_scan_analysis.m
+++ b/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
@@ -1,6 +1,6 @@
kN=2.22;
figure
-ERRBAR = 1; LOGSCALE = 0;
+ERRBAR = 0; LOGSCALE = 0; AU = 1;
resstr={};
msz = 10; lwt = 2.0;
% CO = 'DGGK'; mrkstyl='d';
@@ -66,7 +66,7 @@ for j = 1:numel(directories)
data = {};
for i = 1:N
subdir = subdirectories{i};
- data = compile_results_low_mem(data,subdir,00,10);
+ data = compile_results_low_mem(data,subdir,00,20);
try
Trange = data.Ts0D(end)*[0.5 1.0];
catch % if data does not exist put 0 everywhere
@@ -79,8 +79,8 @@ for j = 1:numel(directories)
data.inputs.K_N = kN;
data.inputs.NU = nus(i);
end
- Trange = data.Ts0D(end)*[0.5 1.0];
- % Trange = [200 400];
+ % Trange = data.Ts0D(end)*[0.5 1.0];
+ Trange = [200 400];
%
[~,it0] = min(abs(Trange(1) -data.Ts0D));
[~,it1] = min(abs(Trange(end)-data.Ts0D));
@@ -93,6 +93,9 @@ for j = 1:numel(directories)
subplot(N,2,2*i-1)
hold on;
Qx = data.HFLUX_X;
+ if AU
+ Qx = Qx./max(Qx);
+ end
T = data.Ts0D;
% Plot heatflux vs time
plot(T,Qx,'DisplayName',[scanvarname,'=',num2str(x(i))],...
diff --git a/wk/p2_heatflux_salpha_convergence.m b/wk/paper_2_scripts_and_results/p2_heatflux_salpha_convergence.m
similarity index 100%
rename from wk/p2_heatflux_salpha_convergence.m
rename to wk/paper_2_scripts_and_results/p2_heatflux_salpha_convergence.m
diff --git a/wk/scan_lin_ITG.m b/wk/paper_2_scripts_and_results/scan_lin_ITG.m
similarity index 100%
rename from wk/scan_lin_ITG.m
rename to wk/paper_2_scripts_and_results/scan_lin_ITG.m
diff --git a/wk/tmp_script.m b/wk/paper_2_scripts_and_results/tmp_script.m
similarity index 100%
rename from wk/tmp_script.m
rename to wk/paper_2_scripts_and_results/tmp_script.m