diff --git a/matlab/load_params.m b/matlab/load_params.m
index c74d7b123c1a9001a71483b193589f82f6ad4787..28ce1ff9cbddfbc31fddb3bb38361c4f18dd0252 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -13,8 +13,8 @@ NZ = h5readatt(filename,'/data/input','nz');
L = h5readatt(filename,'/data/input','Lr');
CLOS = h5readatt(filename,'/data/input','CLOS');
DT_SIM = h5readatt(filename,'/data/input','dt');
-% MU = h5readatt(filename,'/data/input','mu');
-MU = str2num(filename(end-18:end-14)); %bad...
+MU = h5readatt(filename,'/data/input','mu');
+% MU = str2num(filename(end-18:end-14)); %bad...
W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index 8caf105d24c4fc6b9a529f47ae167e7b8812e1da..d4fde7849f639a0695d12bfa3e79e63750db61cf 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -36,7 +36,7 @@ fprintf(fid,[...
'#SBATCH --account=FUA35_TSVVT421\n',...
'#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',...
'module load autoload hdf5 fftw\n',...
-'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KR)]);
+'srun --cpu-bind=cores ./../../../bin/',EXECNAME,' ',num2str(NP_P),' ',num2str(NP_KX)]);
fclose(fid);
system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 746f0e46c5c97fc0fa0a7093ec850c9c67961ed1..cd3d9209c2ff476061112751c19658c0c671d0c1 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -4,10 +4,10 @@ for i_ = 1
%% Load results
if 1% Local results
outfile ='';
-outfile ='HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02';
-% outfile ='HD_study/100x50_L_50_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02';
-% outfile ='kobayashi/100x50_L_50_P_2_J_1_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00';
-% outfile ='v2.7_P_2_J_1/100x50_L_200_P_2_J_1_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_0e+00';
+outfile ='';
+outfile ='';
+outfile ='HD_study/200x100_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02';
+% outfile ='HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-02';
BASIC.RESDIR = ['../results/',outfile,'/'];
BASIC.MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
@@ -19,7 +19,9 @@ outfile ='';
outfile ='';
outfile ='';
outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/out.txt';
+outfile ='';
+outfile ='';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt';
% outfile = outcl{i_};
% load_marconi(outfile);
BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
@@ -87,7 +89,7 @@ for it = 1:numel(Ts2D)
dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
E_turb(it) = sum(sum((1+KR.^2+KZ.^2).*abs(PHI(:,:,it)).^2))- sum((1+kr.^2).*abs(PHI(:,1,it)).^2);
- E_ZF(it) = kr(ikZF)^2*abs(PHI(ikZF,1,it)).^2;
+ E_ZF(it) = sum(kr.^2.*abs(PHI(:,1,it)).^2);
if(W_DENS && W_TEMP)
DENS_E_ = DENS_E(:,:,it); DENS_I_ = DENS_I(:,:,it);
TEMP_E_ = TEMP_E(:,:,it); TEMP_I_ = TEMP_I(:,:,it);
@@ -348,11 +350,11 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
grid on; ylabel('Shear amp.');set(gca,'xticklabel',[]);% legend('show');
subplot(312)
yyaxis left
- plot(Ts2D,SCALE*E_ZF);
- ylabel('ZF energy');
+ plot(Ts2D,SCALE*E_ZF);xlim([0,Ts0D(end)]);
+ ylabel('ZF energy'); ylim([0;1.2*max(SCALE*E_ZF(floor(0.5*numel(Ts2D)):end))]);
yyaxis right
- plot(Ts2D,SCALE*E_turb);
- ylabel('Turb. energy'); ylim([0;1.2*max(SCALE*E_ZF(floor(0.8*numel(Ts2D)):end))]);
+ plot(Ts2D,SCALE*E_turb);xlim([0,Ts0D(end)]);
+ ylabel('Turb. energy'); ylim([0;1.2*max(SCALE*E_turb(floor(0.5*numel(Ts2D)):end))]);
subplot(313)
[TY,TX] = meshgrid(r,Ts2D);
% pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); legend('$\langle \partial_r\phi\rangle_z$') %colorbar;
@@ -361,7 +363,7 @@ fig = figure; FIGNAME = ['ZF_transport_drphi','_',PARAMS];set(gcf, 'Position',
save_figure
end
-if 1
+if 0
%% Space time diagramms
tstart = 0; tend = Ts2D(end);
[~,itstart] = min(abs(Ts2D-tstart));
@@ -370,32 +372,20 @@ trange = itstart:itend;
[TY,TX] = meshgrid(r,Ts2D(trange));
fig = figure; FIGNAME = ['space_time','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(211)
-% pclr = pcolor(TX,TY,squeeze(mean(dens_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
pclr = pcolor(TX,TY,squeeze(mean(ni00(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
shading interp
colormap hot;
caxis([0.0,0.05*max(max(mean(ni00(:,:,its2D:ite2D).*dzphi(:,:,its2D:ite2D),2)))]);
caxis([0.0,0.05]); c = colorbar; c.Label.String ='\langle n_i\partial_z\phi\rangle_z';
xticks([]); ylabel('$r/\rho_s$')
-% legend('Radial part. transport $\langle n_i\partial_z\phi\rangle_z$')
title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
' $\mu_{hd}=$',num2str(MU)]);
-% subplot(312)
-% pclr = pcolor(TX,TY,squeeze(mean(temp_i(:,:,trange).*dzphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
-% shading interp
-% % colormap(bluewhitered(256));
-% xticks([]); ylabel('$r/\rho_s$')
-% legend('Radial part. transport $\langle T_i\partial_z\phi\rangle_z$')
-% title(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
-% ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-% ' $\mu_{hd}=$',num2str(MU)]);
subplot(212)
pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,trange),2))'); set(pclr, 'edgecolor','none'); colorbar;
fieldmax = max(max(mean(abs(drphi(:,:,its2D:ite2D)),2)));
caxis([-fieldmax,fieldmax]); c = colorbar; c.Label.String ='\langle \partial_r\phi\rangle_z';
xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
-% legend('Zonal flow $\langle \partial_r\phi\rangle_z$')
save_figure
end
@@ -416,7 +406,7 @@ end
if 1
%% |phi_k|^2 spectra (Kobayashi 2015 fig 3)
-tstart = 4000; tend = 4000;
+tstart = 2000; tend = 3000;
[~,itstart] = min(abs(Ts2D-tstart));
[~,itend] = min(abs(Ts2D-tend));
trange = itstart:itend;
@@ -439,7 +429,7 @@ ne00_kp = mean(Z_rth,2);
%for theorical trends
a1 = phi_kp(2)*kp_ip(2).^(13/3);
a2 = phi_kp(2)*kp_ip(2).^(3)./(1+kp_ip(2).^2).^(-2);
-fig = figure; FIGNAME = ['cascade','_',PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
+fig = figure; FIGNAME = ['cascade_',num2str(tstart),'to',num2str(tend),PARAMS];set(gcf, 'Position', [100, 100, 500, 500])
% scatter(kperp,phi_k_2,'.k','MarkerEdgeAlpha',0.4,'DisplayName','$|\phi_k|^2$'); hold on; grid on;
plot(kp_ip,phi_kp,'^','DisplayName','$\langle|\phi_k|^2\rangle_{k_\perp}$'); hold on;
plot(kp_ip,ni00_kp, '^','DisplayName','$\langle|N_{i,k}^{00}|^2\rangle_{k_\perp}$'); hold on;
@@ -450,14 +440,14 @@ set(gca, 'XScale', 'log');set(gca, 'YScale', 'log'); grid on
xlabel('$k_\perp \rho_s$'); legend('show','Location','southwest')
title({['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
', $L=',num2str(L),'$, $N=',num2str(Nr),'$'];[' $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$,',...
-' $\mu_{hd}=$',num2str(MU)]});
+' $\mu_{hd}=$',num2str(MU),', $',num2str(tstart),'<t<',num2str(tend),'$']});
xlim([0.1,10]);
save_figure
clear Z_rth phi_kp ni_kp Ti_kp
end
%%
-t0 =000;
+t0 =3000;
[~, it02D] = min(abs(Ts2D-t0));
[~, it05D] = min(abs(Ts5D-t0));
skip_ = 1;
@@ -502,8 +492,8 @@ if 0
GIFNAME = ['ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$n_i^{00}$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
-% create_gif
-create_mov
+create_gif
+% create_mov
end
if 0
%% GC Density electrons
diff --git a/wk/continue_multiple_runs_marconi.m b/wk/continue_multiple_runs_marconi.m
index 7e05b31d7ed287465cfab8b9005e391bab63ee3d..cc2b930440d80c19b0cf337cd5fac5652e318d14 100644
--- a/wk/continue_multiple_runs_marconi.m
+++ b/wk/continue_multiple_runs_marconi.m
@@ -1,5 +1,5 @@
%% Paste the list of continue_run calls
-continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt')
+continue_run('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt')
%% Functions to modify preexisting fort.90 input file and launch on marconi
function [] = continue_run(outfilename)
@@ -44,33 +44,35 @@ function [] = continue_run(outfilename)
line = A{35};
line = line(end-2:end);
if(line(1) == '='); line = line(end); end;
- J2L = str2num(line);
+ J2L = str2num(line)+1;
end
% Change job 2 load in fort.90
- A{35} = [' job2load = ',num2str(J2L+1)];
+ A{35} = [' job2load = ',num2str(J2L)];
disp(A{35})
% Change time step
- A{3} = [' dt = 0.005'];
+ line_= A{3}; dt_old = str2num(line_(12:end));
+ A{3} = [' dt = ',num2str(dt_old)];
% Increase endtime
- A{4} = [' tmax = 10000'];
+ A{4} = [' tmax = 20000'];
% Put non linear term back
- A{41} = [' NL_CLOS = -1'];
+ line_= A{41}; NL_old = str2num(line_(13:end));
+ A{41} = [' NL_CLOS = ',num2str(NL_old)];
% change HD
- line_= A{43};
- mu_old = str2num(line_(13:end));
+ line_= A{43}; mu_old = str2num(line_(13:end));
A{43} = [' mu = ',num2str(mu_old)];
+ % change N
+ line_= A{13}; N_old = str2num(line_(10:end));
+ A{13} = [' Nr = ',num2str(N_old)];
+ A{15} = [' Nz = ',num2str(N_old)];
% change L
- line_= A{14};
- L_old = str2num(line_(8:end));
- A{14} = [' Lx = ',num2str(L_old)];
- A{16} = [' Ly = ',num2str(L_old)];
+ line_= A{14}; L_old = str2num(line_(8:end));
+ A{14} = [' Lr = ',num2str(L_old)];
+ A{16} = [' Lz = ',num2str(L_old)];
% change eta N
- line_= A{53};
- etan_old = str2num(line_(13:end));
+ line_= A{53}; etan_old = str2num(line_(13:end));
A{53} = [' eta_n = ',num2str(etan_old)];
% change eta B
- line_= A{55};
- etab_old = str2num(line_(13:end));
+ line_= A{55}; etab_old = str2num(line_(13:end));
A{55} = [' eta_B = ',num2str(etab_old)];
% Rewrite fort.90
fid = fopen('fort.90', 'w');
diff --git a/wk/load_multiple_outputs_marconi.m b/wk/load_multiple_outputs_marconi.m
index ee1d402194130c2234c24706a2f133d1bb19f287..fc48c1e5bea030b81167d92d47461e9cb9a24028 100644
--- a/wk/load_multiple_outputs_marconi.m
+++ b/wk/load_multiple_outputs_marconi.m
@@ -1,5 +1,7 @@
addpath(genpath('../matlab')) % ... add
%% Paste the list of simulation results to load
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_60_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_1e-03/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_1e-02/out.txt')
-load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e+00_SGGK_CLOS_0_mu_1e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.8_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-03_SGGK_CLOS_0_mu_2e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_1e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/kobayashi/300x150_L_100_P_6_J_3_eta_0.71429_nu_1e-02_PAGK_CLOS_0_mu_0e+00/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_6_J_3/200x100_L_120_P_6_J_3_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt')
+load_marconi('/marconi_scratch/userexternal/ahoffman/HeLaZ/results/v2.7_P_10_J_5/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-02_SGGK_CLOS_0_mu_1e-02/out.txt')
diff --git a/wk/local_run.m b/wk/local_run.m
index 355e9d7c438f7b0f87771367566105a503834f5b..0e002fc07e52a0b70c4a5e7f1b00da5b93a4ec0a 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -7,9 +7,9 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
NU = 0.1; % Collision frequency
ETAB = 0.6; % Magnetic gradient
ETAN = 1.0; % Density gradient
-NU_HYP = 1.0;
+NU_HYP = 0.5;
%% GRID PARAMETERS
-N = 150; % Frequency gridpoints (Nkr = N/2)
+N = 300; % Frequency gridpoints (Nkr = N/2)
L = 100; % Size of the squared frequency domain
P = 2;
J = 1;
@@ -17,7 +17,7 @@ MU_P = 0.0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0.0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
TMAX = 1000; % Maximal time unit
-DT = 1e-2; % Time step
+DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS5D = 1/200; % Sampling per time unit for 5D arrays
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 1e8aabdb1f46c1265a404c20ea07b1a02fefa3d0..690b7f670a2c1e548e5ae5963662d1d030d41810 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,8 +1,8 @@
clear all;
addpath(genpath('../matlab')) % ... add
-SUBMIT = 1; % To submit the job automatically
+SUBMIT = 0; % To submit the job automatically
% EXECNAME = 'helaz_dbg';
- EXECNAME = 'helaz_2.73';
+ EXECNAME = 'helaz_2.8';
for ETAB = [0.6]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
@@ -15,7 +15,7 @@ if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
CLUSTER.MEM = '128GB'; % Memory
CLUSTER.JNAME = 'HeLaZ';% Job name
NP_P = 2; % MPI processes along p
-NP_KR = 24; % MPI processes along kr
+NP_KX = 24; % MPI processes along kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 1e-3; % Collision frequency
@@ -81,7 +81,7 @@ ETAT = 0.0; % Temperature gradient
ETAN = 1.0; % Density gradient
TAU = 1.0; % e/i temperature ratio
% Compute processes distribution
-Ntot = NP_P * NP_KR;
+Ntot = NP_P * NP_KX;
Nnodes = ceil(Ntot/48);
Nppn = Ntot/Nnodes;
CLUSTER.NODES = num2str(Nnodes); % MPI process along p
@@ -91,7 +91,7 @@ CLUSTER.CPUPT = '1'; % CPU per task
setup
write_sbash_marconi
system('rm fort.90 setup_and_run.sh batch_script.sh');
-if(mod(NP_P*NP_KR,48)~= 0)
+if(mod(NP_P*NP_KX,48)~= 0)
disp('WARNING : unused cores (ntot cores must be a 48 multiple)');
end
if(SUBMIT)
diff --git a/wk/photomaton.m b/wk/photomaton.m
index 32b5454dbc624c21005b9ca6d71f1714758b52e2..b8801d6ef1ebbff96a0f7890dbb4e45a308b9432 100644
--- a/wk/photomaton.m
+++ b/wk/photomaton.m
@@ -1,19 +1,19 @@
if 0
%% Photomaton : real space
% FIELD = ni00; FNAME = 'ni00'; FIELDLTX = '$n_i^{00}$'; XX = RR; YY = ZZ;
-FIELD = ne00; FNAME = 'ne00'; FIELDLTX = '$n_e^{00}$'; XX = RR; YY = ZZ;
+% FIELD = ne00; FNAME = 'ne00'; FIELDLTX = '$n_e^{00}$'; XX = RR; YY = ZZ;
% FIELD = dens_i; FNAME = 'ni'; FIELDLTX = '$n_i$'; XX = RR; YY = ZZ;
% FIELD = dens_e; FNAME = 'ne'; FIELDLTX = '$n_e$'; XX = RR; YY = ZZ;
% FIELD = temp_i; FNAME = 'Ti'; FIELDLTX = '$T_i$'; XX = RR; YY = ZZ;
% FIELD = temp_e; FNAME = 'Te'; FIELDLTX = '$T_e$'; XX = RR; YY = ZZ;
-% FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
+FIELD = phi; FNAME = 'phi'; FIELDLTX = '$\phi$'; XX = RR; YY = ZZ;
% FIELD = drphi; FNAME = 'ZF'; FIELDLTX = '$u^{ZF}_z$'; XX = RR; YY = ZZ;
% FIELD = -dr2phi; FNAME = 'shear'; FIELDLTX = '$s$'; XX = RR; YY = ZZ;
TNAME = [];
-tf = 250; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 450; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 1200; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
-tf = 1500; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 2500; [~,it1] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 3000; [~,it2] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 3500; [~,it3] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
+tf = 4000; [~,it4] = min(abs(Ts2D-tf)); TNAME = [TNAME,'_',num2str(tf)];
fig = figure; FIGNAME = [FNAME,TNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 350])
plt = @(x) x./max(max(x));
subplot(141)