diff --git a/matlab/assemble_cosolver_matrices.m b/matlab/assemble_cosolver_matrices.m
index 6c0afb15853c1cde9e9410f546d2b905709035a6..35c03c69c8a45040a68e97c082a3ee4452caf7cf 100644
--- a/matlab/assemble_cosolver_matrices.m
+++ b/matlab/assemble_cosolver_matrices.m
@@ -1,6 +1,6 @@
codir = '/home/ahoffman/cosolver/';
-matname= 'gk_landauii_P16_J9_dk_5e-2_km_2.0_NFLR_8';
-matdir = [codir,'LDGKii_P16_J9_NFLR_8/matrices/'];
+matname= 'gk_landau_P16_J9_dk_5e-2_km_2.0_NFLR_8';
+matdir = [codir,'LDGK_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 = ['self.',n_string,'.h5'];
+ filename = ['ei.',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 = '/Caapj/Ciipj'; infile = [matdir,'self.',n_string,'.h5'];
-% Pmaxe = h5readatt(infile,olddname,'Pmaxe');
-% Jmaxe = h5readatt(infile,olddname,'Jmaxe');
+olddname = '/Ceipj/CeipjF'; infile = [matdir,'ei.',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);
%
@@ -68,4 +68,4 @@ h5create(outfile,'/coordkperp',numel(kperp(1:Nmax)));
h5write (outfile,'/coordkperp',kperp(1:Nmax)');
fid = H5F.open(outfile);
-H5F.close(fid);
\ No newline at end of file
+H5F.close(fid);
diff --git a/matlab/compute/ballooning_representation.m b/matlab/compute/ballooning_representation.m
new file mode 100644
index 0000000000000000000000000000000000000000..dbac83f8f8698895f66e5a272434c9fcee460207
--- /dev/null
+++ b/matlab/compute/ballooning_representation.m
@@ -0,0 +1,46 @@
+function [field_chi,chi] = ballooning_representation(field,shear,Npol,kx,iky,ky,z)
+ dims = size(field);
+ Nkx = dims(2);
+ % Apply ballooning transform
+ if(Nkx > 1)
+ nexc = round(ky(2)*shear*2*pi/kx(2));
+ else
+ nexc = 1;
+ end
+ is = max(1,iky-1);
+ Npi = (Nkx-1)-2*nexc*(is-1);
+ if(Npi <= 1)
+ ordered_ikx = 1;
+ elseif(mod(Nkx,2) == 0)
+ tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
+ ordered_ikx = [tmp_(end:-1:1), 1:is:(Nkx/2)];
+ else
+ Np_ = (Nkx+1)/(2*is);
+ ordered_ikx = [(Np_+1):Nkx 1:Np_];
+ end
+ try
+ idx=kx./kx(2);
+ catch
+ idx=0;
+ end
+ idx= idx(ordered_ikx);
+ Nkx = numel(idx);
+
+ field_chi = zeros( Nkx*dims(3) ,1);
+ chi = zeros( Nkx*dims(3) ,1);
+
+ for i_ =1:Nkx
+ start_ = (i_-1)*dims(3) +1;
+ end_ = i_*dims(3);
+ ikx = ordered_ikx(i_);
+ field_chi(start_:end_) = field(iky,ikx,:);
+ end
+
+ % Define ballooning angle
+ for i_ =1: Nkx
+ for iz=1:dims(3)
+ ii = dims(3)*(i_-1) + iz;
+ chi(ii) =z(iz) + 2*pi*Npol*idx(i_)./is;
+ end
+ end
+end
diff --git a/matlab/compute/invert_kxky_to_kykx_gene_results.m b/matlab/compute/invert_kxky_to_kykx_gene_results.m
index 3f320d653472aa52e6583aee1097355680323f1e..69d3cbb07d4d3592630722c46cd022ac7688f5b1 100644
--- a/matlab/compute/invert_kxky_to_kykx_gene_results.m
+++ b/matlab/compute/invert_kxky_to_kykx_gene_results.m
@@ -17,8 +17,11 @@ if numel(data.grids.kx)>1
else
dkx = 1;
end
-
-dky = data.grids.ky(2);
+if data.grids.Nky > 1
+ dky = data.grids.ky(2);
+else
+ dky = data.grids.ky(1);
+end
Lx = 2*pi/dkx; Ly = 2*pi/dky;
x = linspace(-Lx/2,Lx/2,data.grids.Nx+1); x = x(1:end-1);
y = linspace(-Ly/2,Ly/2,data.grids.Ny+1); y = y(1:end-1);
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 907aa85e2122328d82dcee05e23763fb089c5457..921dd74352eea0f4de251a2a142816e5ada44574 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -83,6 +83,7 @@ for i = 1:2
amp = MODES;
i_=1;
for ik = MODES
+
to_measure = log(plt(field,ik));
gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
gamma(i_) = gr(1);
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
index 196e899617f3ed7af330823bf805d1bb3f7e4c1a..7ebeabcd2a72769c6942e7851c077609a4b61355 100644
--- a/matlab/load/load_gene_data.m
+++ b/matlab/load/load_gene_data.m
@@ -1,6 +1,10 @@
function [ DATA ] = load_gene_data( folder )
%to load gene data as for HeLaZ results
-namelist = read_namelist([folder,'parameters']);
+try
+ namelist = read_namelist([folder,'parameters']);
+catch
+ namelist = read_namelist([folder,'parameters.dat']);
+end
DATA.namelist = namelist;
DATA.folder = folder;
%% Grid
@@ -29,7 +33,11 @@ if numel(DATA.grids.kx)>1
else
dkx = 1;
end
-dky = DATA.grids.ky(2);
+if DATA.grids.Nky > 1
+ dky = DATA.grids.ky(2);
+else
+ dky = DATA.grids.ky(1);
+end
Lx = 2*pi/dkx; Ly = 2*pi/dky;
x = linspace(-Lx/2,Lx/2,DATA.grids.Nx+1); x = x(1:end-1);
y = linspace(-Ly/2,Ly/2,DATA.grids.Ny+1); y = y(1:end-1);
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 191cdf1a6f2650d4a6cc4f5c5437b4f7d9e298d4..ce0f179434b154b96307437fcf4d9731898bce5c 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -37,13 +37,20 @@ 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_I = h5readatt(filename,'/data/input/model','tau_i');
+ try
+ DATA.inputs.TAU_I = h5readatt(filename,'/data/input/model','tau_i');
+ catch
+ DATA.inputs.TAU_I = 1.0;
+ end
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';
-
+ try
+ DATA.inputs.ADIAB_I = h5readatt(filename,'/data/input/model','ADIAB_I') == 'y';
+ catch
+ DATA.inputs.ADIAB_I = 0;
+ end
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';
DATA.inputs.W_NA00 = h5readatt(filename,'/data/input/diag_par','write_Na00') == 'y';
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 0b8eaea271fb1efa6af82169e6a62a7dc47ddab5..299c9b1119d5a8d2a7b5393db32b249fb05a2455 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -4,6 +4,7 @@ function [FIG] = plot_ballooning(data,options)
[~,it1] = min(abs(data.Ts3D - options.time_2_plot(end)));
[~,ikyarray] = min(abs(data.grids.ky - options.kymodes));
ikyarray = unique(ikyarray);
+ dz = data.grids.z(2) - data.grids.z(1);
phi_real=real(data.PHI(:,:,:,it1));
phi_imag=imag(data.PHI(:,:,:,it1));
ncol = 1;
@@ -12,55 +13,19 @@ function [FIG] = plot_ballooning(data,options)
psi_imag=imag(data.PSI(:,:,:,it1));
ncol = 2;
end
- % Apply ballooning transform
- if(data.grids.Nkx > 1)
- nexc = round(data.grids.ky(2)*data.inputs.SHEAR*2*pi/data.grids.kx(2));
- else
- nexc = 1;
+ if options.PLOT_KP
+ ncol = 3;
end
+ % % Apply ballooning transform
+ % if(data.grids.Nkx > 1)
+ % nexc = round(data.grids.ky(2)*data.inputs.SHEAR*2*pi/data.grids.kx(2));
+ % else
+ % nexc = 1;
+ % end
for iky=ikyarray
- dims = size(phi_real);
- Nkx = dims(2);
- is = max(1,iky-1);
- Npi = (Nkx-1)-2*nexc*(is-1);
- if(Npi <= 1)
- ordered_ikx = 1;
- elseif(mod(Nkx,2) == 0)
- tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
- ordered_ikx = [tmp_(end:-1:1), 1:is:(Nkx/2)];
- else
- Np_ = (Nkx+1)/(2*is);
- ordered_ikx = [(Np_+1):Nkx 1:Np_];
- end
- try
- idx=data.grids.kx./data.grids.kx(2);
- catch
- idx=0;
- end
- idx= idx(ordered_ikx);
- Nkx = numel(idx);
-
- phib_real = zeros( Nkx*dims(3) ,1);
- phib_imag = phib_real;
- b_angle = phib_real;
-
- for i_ =1:Nkx
- start_ = (i_-1)*dims(3) +1;
- end_ = i_*dims(3);
- ikx = ordered_ikx(i_);
- phib_real(start_:end_) = phi_real(iky,ikx,:);
- phib_imag(start_:end_) = phi_imag(iky,ikx,:);
- end
- % Define ballooning angle
- 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*data.grids.Npol*idx(i_)./is;
- end
- end
-
+ [phib_real,b_angle] = ballooning_representation(phi_real,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ [phib_imag,~] = ballooning_representation(phi_imag,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
phib = phib_real(:) + 1i * phib_imag(:);
% normalize real and imaginary parts at chi =0
if options.normalized
@@ -82,19 +47,17 @@ function [FIG] = plot_ballooning(data,options)
ylabel('$\phi_B (\chi)$');
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)),']$']);
- for i_ = 2*[1:data.grids.Nkx]-(data.grids.Nkx+1)
- xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ % z domain start end point
+ for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ end
+ for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_)-dz/pi,'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;
- end_ = i_*dims(3);
- ikx = ordered_ikx(i_);
- psib_real(start_:end_) = psi_real(iky,ikx,:);
- psib_imag(start_:end_) = psi_imag(iky,ikx,:);
- end
+ [psib_real,b_angle] = ballooning_representation(psi_real,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ [psib_imag,~] = ballooning_representation(psi_imag,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+
psib = psib_real(:) + 1i * psib_imag(:);
psib_norm = psib / normalization;
psib_real_norm = real( psib_norm);
@@ -104,15 +67,37 @@ 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)
+ % z domain start end point
+ for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
end
+ for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ end
xlabel('$\chi / \pi$')
ylabel('$\psi_B (\chi)$');
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
+ if options.PLOT_KP
+ kperp = h5read(data.outfilenames{1},'/data/grid/coordkp');
+ [kperpb,b_angle] = ballooning_representation(kperp,data.inputs.SHEAR,data.grids.Npol,data.grids.kx,iky,data.grids.ky,data.grids.z);
+ subplot(numel(ikyarray),ncol,ncol*(iplot-1)+3)
+ plot(b_angle/pi, kperpb,'-k'); hold on;
+ % z domain start end point
+ for i_ = 2*[1:data.grids.Nkx-1]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_),'HandleVisibility','off'); hold on;
+ end
+ for i_ = 2*[2:data.grids.Nkx]-(data.grids.Nkx+1)
+ xline(data.grids.Npol*(i_)-dz/pi,'HandleVisibility','off'); hold on;
+ end
+ xlabel('$\chi / \pi$')
+ ylabel('$k_\perp (\chi)$');
+ 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
iplot = iplot + 1;
end
end
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index 8a813cff4501b6cd5cff46fccc8874a82b8e7e0a..a9afcf846fa39a7f3b64c00a69a0512d2442919d 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -98,38 +98,25 @@ if 0
%%
mfns = {...
'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';...
-% '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5';...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5';...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5';...
-% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5';...
+ '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';...
'/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';...
-% '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5';...
'/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_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_landau_P16_J9_dk_5e-2_km_2.0_NFLR_8.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';...
};
CONAME_A = {...
'SG 20 10';...
-% 'PA 20 10';...
-% 'FC 10 5 NFLR 4';...
-% 'FC 10 5 NFLR 12';...
-% 'FC 10 5 NFLR 12 k<2';, ...
+ 'PA 20 10';...
'LD 6 3 NFLR 20'; ...
-% 'FC 4 2 NFLR 6';...
'FC 4 2 NFLR 12'; ...
'LD 10 5 NFLR 12'; ...
'LD 11 7 NFLR 16'; ...
'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'; ...
+ 'LD 16 9 NFLR 8'; ...
% 'Hacked SG A';...
% 'Hacked SG B';...
};
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index d1caeced8f68f96b0da031403ae1ce41fe642f3d..5fb56042ae281f23a85fa54d7cc09f79419810fd 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -13,7 +13,7 @@ for i_ = 1:numel(names)
end
NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
if NPLOT > 0
- fig = figure;
+ fig = figure;
if options.SHOW_METRICS
subplot(3,NPLOT,1*NPLOT)
for i = 1:6
@@ -32,6 +32,7 @@ if NPLOT > 0
plot(data.grids.z, geo_arrays(:,i),'DisplayName',names{i}); hold on;
end
xlim([min(data.grids.z),max(data.grids.z)]); legend('show');
+ legend('Interpreter','none')
end
if options.SHOW_FLUXSURF
subplot(1,2,1)
diff --git a/matlab/plot/show_moments_spectrum.m b/matlab/plot/show_moments_spectrum.m
index a33afa4e2e294e8f2fea3b762254b4a86dd69fd2..ab9d3aea55bd247e410291c6560aa38f823805d4 100644
--- a/matlab/plot/show_moments_spectrum.m
+++ b/matlab/plot/show_moments_spectrum.m
@@ -11,10 +11,10 @@ if OPTIONS.ST == 0
end
if OPTIONS.LOGSCALE
logname = 'log';
- compress = @(x,ia) log(sum(abs((x(:,:,:,:,:))),4));
+ compress = @(x,ia) log(sum(abs((x(ia,:,:,:,:))),4));
else
logname = '';
- compress = @(x,ia) sum(abs((x(:,:,:,:,:))),4);
+ compress = @(x,ia) sum(abs((x(ia,:,:,:,:))),4);
end
for ia = 1:DATA.inputs.Na
Napjz = compress(DATA.Napjz,ia);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 555931e62f8ddb7193996542ce716d03ac91c811..e50a65d18b4e58a20720ecc26c830bf8713d5944 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -29,11 +29,11 @@ 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,[' 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,[' 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,[' 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']);
@@ -87,7 +87,7 @@ fprintf(fid,'/\n');
if(strcmp(MODEL.ADIAB_I,'.false.'))
fprintf(fid,'&SPECIES\n');
- fprintf(fid, ' name_ = ions \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']);
@@ -97,7 +97,7 @@ if(strcmp(MODEL.ADIAB_I,'.false.'))
end
if(strcmp(MODEL.ADIAB_E,'.false.'))
fprintf(fid,'&SPECIES\n');
- fprintf(fid, ' name_ = electrons \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/AUG_ky_PJ_scan.m b/wk/AUG_ky_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..2b562a895e1d8c90af3472024071db00a88bbe74
--- /dev/null
+++ b/wk/AUG_ky_PJ_scan.m
@@ -0,0 +1,224 @@
+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_AUG_scan'; % Name of the simulation
+RERUN = 1; % rerun if the data does not exist
+RUN = 1;
+K_Ne = 34; % ele Density '''
+K_Te = 56; % ele Temperature '''
+K_Ni = 34; % ion Density gradient drive
+K_Ti = 21; % 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 = 1.52e-4; % electron plasma beta
+MHD_PD = 0; % MHD pressure drift
+% Geometry
+GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
+% 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.38; % 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= 'miller';
+ EPS = 0.3; % inverse aspect ratio
+ Q0 = 5.7; % safety factor
+ SHEAR = 4.6;%*EPS; % magnetic shear
+ KAPPA = 1.8; % elongation
+ S_KAPPA = 0.0;
+ DELTA = 0.4; % 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
+ PB_PHASE= 0;
+ MHD_PD = 0; % MHD pressure drift
+ %% 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;
+ ADIAB_I = 0; % adiabatic ion model
+ %% 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,' 1 2 2 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_ky_PJ_scan.m b/wk/KBM_ky_PJ_scan.m
index 1fe7d44587a1e1b23b7d3e3b6cc7f25e8dca716e..45aff676d44927d3799d9689154b55756a104b68 100644
--- a/wk/KBM_ky_PJ_scan.m
+++ b/wk/KBM_ky_PJ_scan.m
@@ -10,7 +10,7 @@ EXECNAME = 'gyacomo23_sp';
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
+RERUN = 0; % rerun if the data does not exist
RUN = 1;
K_Ne = 3; % ele Density '''
K_Te = 4.5; % ele Temperature '''
@@ -29,8 +29,8 @@ NU = 1e-2;
NA = 2;
BETA = 0.03; % electron plasma beta
% Geometry
-% GEOMETRY= 'miller';
-GEOMETRY= 's-alpha';
+GEOMETRY= 'miller';
+% GEOMETRY= 's-alpha';
SHEAR = 0.8; % magnetic shear
% time and numerical grid
DT0 = 5e-3;
diff --git a/wk/KBM_load_metadata_scan.m b/wk/KBM_load_metadata_scan.m
index 5309527b6eff7ba71007b7e46860e9222313c21f..b86161b30dce5a2e4937409ecc781439b1afbecc 100644
--- a/wk/KBM_load_metadata_scan.m
+++ b/wk/KBM_load_metadata_scan.m
@@ -6,8 +6,9 @@ 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';
-
+% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
+% datafname = 'lin_KBM/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.03.mat';
+datafname ='lin_AUG_scan/12x24_ky_0.05_0.75_P_2_16_DGGK_0.01_be_0.000152.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 0;
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 20777d70c099eddf7dd18ec7338862f53b312593..15c758bd6e471fc8ac9656c75125921f23067977 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -62,13 +62,15 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/kT_scan_nu0/30x16x128x64x24/kT_4.0/';
% folder = '/misc/gene_results/kT_scan_nu0/42x24x128x64x24/kT_6.5/';
% Miller CBC
-folder = '/misc/gene_results/Miller_KT_6.96_128x64x24x16x8/output/';
-% folder = '/misc/gene_results/CBC/Miller_KT_4.1_128x64x24x16x8/';
+% folder = '/misc/gene_results/Miller_AUG/Npol_1/output/';
+% folder = '/misc/gene_results/Miller_AUG/Npol_5/output/';
+% folder = '/misc/gene_results/Miller_AUG/circ/';
+folder = '/misc/gene_results/Miller_AUG/linear/Npol_1/';
% debug ? shearless
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_00/';
% folder = '/misc/gene_results/CBC/shearless_CBC_128x64x24x24x12_01/';
-if 0
+if 1
%% FULL DATA LOAD (LONG)
gene_data = load_gene_data(folder);
gene_data.FIGDIR = folder;
@@ -85,7 +87,7 @@ if 0
dashboard(gene_data);
end
-if 1
+if 0
%% ONLY HEAT FLUX
nrgfile = 'nrg.dat.h5';
% nrgfile = 'nrg_1.h5';
@@ -168,7 +170,7 @@ gene_data.a = data.EPS * 2000;
create_film(gene_data,options,'.gif')
end
-if 0
+if 1
%% Geometry
names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
'$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
diff --git a/wk/fast_analysis.m b/wk/fast_analysis.m
index 0c58560405fa2baf611637c4d187dd78a91cf186..7e73e23389ea84ee0b4915b3b79259f9ff17468d 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';
@@ -90,6 +90,9 @@ resdir = '/home/ahoffman/gyacomo/results/dev/3D_kine_zpinch_test';
% resdir = '/misc/gyacomo23_outputs/paper_2_GYAC23/kT_scan_nu_1e-3/9x5x128x64x24/kT_7.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';
+
+%% Paper 3
+resdir = '/misc/gyacomo23_outputs/paper_3_GYAC23/nonlin_KBM';
%%
J0 = 00; J1 = 20;
@@ -156,9 +159,10 @@ options.AXISEQUAL = 0;
options.NORMALIZE = 0;
% options.NAME = 'N_i^{00}';
options.NAME = '\phi';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
options.COMP = 'avg';
-options.TIME = [100 200 300];
+options.TIME = [0:0.02:0.14];
+% options.TIME = data.Ts3D(1:2:end);
options.RESOLUTION = 256;
fig = photomaton(data,options);
% save_figure(data,fig)
diff --git a/wk/lin_AUG.m b/wk/lin_AUG.m
index ae21fb232f1cdfe21f58f263c1e5b65a76418c0e..8cfded1e89c2b8bb6f06fbc5abfb222e00f48ff3 100644
--- a/wk/lin_AUG.m
+++ b/wk/lin_AUG.m
@@ -32,38 +32,34 @@ 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
+NPOL = 1; % Number of poloidal turns
%% Set up grid parameters
-P = 2;
+P = 4;
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
+NX = 6; % real space x-gridpoints
+NY = 8; % 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)
+LY = 2*pi/0.4; % Size of the squared frequency domain in y direction
+NZ = 32*NPOL; % 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= 'circular';
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
+SHEAR = 4.6;%*EPS; % 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
@@ -85,7 +81,7 @@ 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
+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
@@ -117,6 +113,7 @@ 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
+ADIAB_I = 0; % adiabatic ion model
%%-------------------------------------------------------------------------
%% RUN
@@ -125,10 +122,10 @@ setup
% 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 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;'];
+ % RUN =['time mpirun -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
end
@@ -169,9 +166,18 @@ 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.kymodes = [0.1:0.1:1];
options.normalized = 1;
+options.PLOT_KP = 0;
% options.field = 'phi';
fig = plot_ballooning(data,options);
+title(GEOMETRY)
end
+if 0
+%% plot geom and metric coefficients
+options.SHOW_FLUXSURF = 1;
+options.SHOW_METRICS = 1;
+[fig, geo_arrays] = plot_metric(data,options);
+title(GEOMETRY)
+end
\ No newline at end of file
diff --git a/wk/lin_ETG_adiab_i.m b/wk/lin_ETG_adiab_i.m
new file mode 100644
index 0000000000000000000000000000000000000000..1c3d907e7082d74f0629cdd4f99507184b42f221
--- /dev/null
+++ b/wk/lin_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.005; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 2.22; % ele Density '''
+K_Te = 6.96; % ele Temperature '''
+K_Ni = 2.22; % ion Density gradient drive
+K_Ti = 6.96; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = 0; % adiabatic electron model
+ADIAB_I = 0; % adiabatic ion model
+BETA = 0.001; % 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 = 8; % real space x-gridpoints
+NY = 2; % real space y-gridpoints
+LX = 2*pi/0.05; % Size of the squared frequency domain in x direction
+LY = 2*pi/10.5; % 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';
+% GEOMETRY= 'z-pinch';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % 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 = 50; % Maximal time unit
+DT = 1e-4; % Time step
+DTSAVE0D = 1; % Sampling per time unit for 0D arrays
+DTSAVE2D = -1; % Sampling per time unit for 2D arrays
+DTSAVE3D = 1; % 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 = '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.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 = 1.0e-8; % Initial noise amplitude
+BCKGD0 = 0.0e-8; % 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 3 2 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/lin_Entropy.m b/wk/lin_Entropy.m
new file mode 100644
index 0000000000000000000000000000000000000000..6197a97f2f5691770d7a97e6eaf2453d2071c865
--- /dev/null
+++ b/wk/lin_Entropy.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.00; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 2.5; % ele Density '''
+K_Te = K_Ne/4; % ele Temperature '''
+K_Ni = 2.5; % ion Density gradient drive
+K_Ti = K_Ni/4; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+NA = 2; % number of kinetic species
+ADIAB_E = 0; % adiabatic electron model
+ADIAB_I = 0; % adiabatic ion model
+BETA = 0.000; % electron plasma beta
+MHD_PD = 0; % MHD pressure drift
+%% Set up grid parameters
+P = 4;
+J = P/2;%P/2;
+PMAX = P; % Hermite basis size
+JMAX = J; % Laguerre basis size
+NX = 2; % real space x-gridpoints
+NY = 12; % real space y-gridpoints
+LX = 2*pi/0.05; % 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';
+GEOMETRY= 'z-pinch';
+EPS = 0.18; % inverse aspect ratio
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % 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 = 100; % 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 = 1; % 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 = '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.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 = 1.0e-4; % Initial noise amplitude
+BCKGD0 = 0.0e-8; % 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 4 1 0;'];
+ % RUN =['time mpirun -np 6 ',gyacomodir,'bin/',EXECNAME,' 1 3 2 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/lin_ITG_salpha.m b/wk/lin_ITG_salpha.m
index 8616c009349c73bd3b1fdb1b4a430321a350b411..26a99fb7638511514802facf8b1dc19bb8b5ee0d 100644
--- a/wk/lin_ITG_salpha.m
+++ b/wk/lin_ITG_salpha.m
@@ -25,7 +25,7 @@ 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 = 5.3; % ion Temperature
+K_Ti = 6.96; % 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
@@ -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 = 1; % Gyrokinetic operator
+GKCO = 0; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
HRCY_CLOS = 'truncation'; % Closure model for higher order moments
diff --git a/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m b/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
index 60627f79cbdf77bd3abe8fd568748618087381d6..a3223ace04e9c2f75821de934d56b96096dbe179 100644
--- a/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
+++ b/wk/paper_2_scripts_and_results/nonlin_nu_scan_analysis.m
@@ -1,11 +1,11 @@
kN=2.22;
figure
-ERRBAR = 0; LOGSCALE = 0; AU = 1;
+ERRBAR = 0; LOGSCALE = 0; AU = 0;
resstr={};
msz = 10; lwt = 2.0;
% CO = 'DGGK'; mrkstyl='d';
-CO = 'SGGK'; mrkstyl='s';
-% CO = 'LDGK'; mrkstyl='o';
+% CO = 'SGGK'; mrkstyl='s';
+CO = 'LDGK'; mrkstyl='o';
% GRAD = 'kT_7.0'; kT = 6.96;
GRAD = 'kT_5.3'; kT = 5.3;
% GRAD = 'kT_4.5'; kT = 4.5;
@@ -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));