From 6e9d3d51b66cb5e9d8f97176c2ccb047ad3d1011 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 20 Apr 2021 10:22:43 +0200
Subject: [PATCH] scripts update
---
matlab/compile_results.m | 14 +++++++
matlab/create_gif_1D_phi.m | 62 +++++++++++++++++++++++++++++++
matlab/profiler.m | 72 ++++++++++++++++++++++++++----------
matlab/write_sbash_daint.m | 3 +-
matlab/write_sbash_marconi.m | 5 +--
wk/analysis_2D.m | 72 +++++++++++++++++++++---------------
wk/compute_collision_mat.m | 12 +++---
wk/daint_run.m | 58 ++++++++++++++++++++---------
wk/local_run.m | 32 ++++++++--------
wk/marconi_run.m | 35 +++++++++---------
10 files changed, 256 insertions(+), 109 deletions(-)
create mode 100644 matlab/create_gif_1D_phi.m
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 3e4bb910..a2689bc1 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -1,5 +1,11 @@
CONTINUE = 1;
JOBNUM = 0; JOBFOUND = 0;
+TJOB_SE = []; % Start and end times of jobs
+NU_EVOL = []; % evolution of parameter nu between jobs
+ETAB_EVOL= []; %
+L_EVOL = []; %
+
+% FIELDS
Nipj_ = []; Nepj_ = [];
Ni00_ = []; Ne00_ = [];
GGAMMA_ = [];
@@ -68,6 +74,13 @@ while(CONTINUE)
Sepj_ = cat(5,Sepj_,Sepj);
end
+ % Evolution of simulation parameters
+ load_params
+ TJOB_SE = [TJOB_SE Ts0D(1) Ts0D(end)];
+ NU_EVOL = [NU_EVOL NU NU];
+ ETAB_EVOL = [ETAB_EVOL ETAB ETAB];
+ L_EVOL = [L_EVOL L L];
+
JOBFOUND = JOBFOUND + 1;
LASTJOB = JOBNUM;
elseif (JOBNUM > 20)
@@ -77,6 +90,7 @@ while(CONTINUE)
JOBNUM = JOBNUM + 1;
Pe_old = Pe_new; Je_old = Je_new;
Pi_old = Pi_new; Ji_old = Ji_new;
+
end
GGAMMA_RI = GGAMMA_; PGAMMA_RI = PGAMMA_; Ts0D = Ts0D_;
Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
diff --git a/matlab/create_gif_1D_phi.m b/matlab/create_gif_1D_phi.m
new file mode 100644
index 00000000..f68761cf
--- /dev/null
+++ b/matlab/create_gif_1D_phi.m
@@ -0,0 +1,62 @@
+title1 = GIFNAME;
+FIGDIR = BASIC.RESDIR;
+if ~exist(FIGDIR, 'dir')
+ mkdir(FIGDIR)
+end
+
+GIFNAME = [FIGDIR, GIFNAME,'.gif'];
+
+% Set colormap boundaries
+hmax = max(max(max(FIELD)));
+hmin = min(min(min(FIELD)));
+fieldabsmin = min(min(abs(FIELD(:,:))));
+fieldabsmax = max(max(abs(FIELD(:,:))));
+flag = 0;
+if hmax == hmin
+ disp('Warning : h = hmin = hmax = const')
+else
+% Setup figure frame
+fig = figure;
+subplot(211)
+ plot(X,FIELD(:,1)); % to set up
+ axis tight manual % this ensures that getframe() returns a consistent size
+ in = 1;
+ nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
+ for n = FRAMES % loop over selected frames
+ scale = max(FIELD(:,n))*SCALING + (1-SCALING);
+
+ subplot(211)
+ plot(X,FIELD(:,n)/scale,linestyle);
+ if (YMIN ~= YMAX && XMIN ~= XMAX)
+ ylim([YMIN,YMAX]); xlim([XMIN,XMAX]);
+ end
+ title(['$t \approx$', sprintf('%.3d',ceil(T(n))), ', scaling = ',sprintf('%.1e',scale)]);
+ xlabel(XNAME); ylabel(FIELDNAME);
+ drawnow
+ % Capture the plot as an image
+ frame = getframe(fig);
+ im = frame2im(frame);
+ [imind,cm] = rgb2ind(im,64);
+ % Write to the GIF File
+ if in == 1
+ imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf);
+ else
+ imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY);
+ end
+ % terminal info
+ while nbytes > 0
+ fprintf('\b')
+ nbytes = nbytes - 1;
+ end
+ nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+ in = in + 1;
+
+ subplot(212)
+ ylim([fieldabsmin,fieldabsmax]); xlim([T(FRAMES(1)),T(FRAMES(end))]); hold on;
+ plot(T(n),abs(max(FIELD(:,n))),'.k'); xlabel('$\rho_s/c_s$')
+
+ end
+ disp(' ')
+ disp(['Gif saved @ : ',GIFNAME])
+end
+
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 37d363ba..57b84701 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -29,14 +29,26 @@ TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/10,10]);
TIME_PER_STEP = sum(TIME_PER_FCT,2);
TIME_PER_CPU = trapz(Ts0D(2:end),TIME_PER_STEP);
+rhs_Ta = mean(diff(rhs_Tc));
+adv_field_Ta = mean(diff(adv_field_Tc));
+ghost_Ta = mean(diff(ghost_Tc));
+clos_Ta = mean(diff(clos_Tc));
+coll_Ta = mean(diff(coll_Tc));
+poisson_Ta = mean(diff(poisson_Tc));
+Sapj_Ta = mean(diff(Sapj_Tc));
+checkfield_Ta = mean(diff(checkfield_Tc));
+diag_Ta = mean(diff(diag_Tc));
+
+NSTEP_PER_SAMP= mean(diff(Ts0D))/DT_SIM;
+
%% Plots
-TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(poisson_Tc); diff(ghost_Tc);...
- diff(Sapj_Tc); diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
-TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/8,8]);
+if 1
+ %% Area plot
fig = figure;
p1 = area(Ts0D(2:end),TIME_PER_FCT,'LineStyle','none');
-legend('Compute RHS','Adv. fields','Poisson', 'comm', 'Sapj','Check+Sym','Diag','Missing')
+for i = 1:10; p1(i).FaceColor = rand(1,3); end;
+legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
xlim([Ts0D(2),Ts0D(end)]);
title(sprintf('Proc. 1, total sim. time ~%.0f [h]',CPUTIME/3600))
@@ -44,18 +56,40 @@ hold on
FIGNAME = 'profiler';
save_figure
-%% Plots
-% fig = figure;
-%
-% p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none');
-% legend('Compute RHS','Adv. fields','Poisson','Sapj','Check+Sym','Diag','Missing')
-% xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
-% ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
-% hold on
-% yyaxis right
-% p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0);
-% ylabel('Step Comp. Time [s]')
-% ylim([0,1.1*max(diff(total_Tc))])
-% set(gca,'ycolor','r')
-% FIGNAME = 'profiler';
-% save_figure
\ No newline at end of file
+else
+ %% Normalized Area plot
+fig = figure;
+
+p1 = area(Ts0D(2:end),100*TIME_PER_FCT./diff(total_Tc),'LineStyle','none', 'FaceColor','flat');
+for i = 1:10; p1(i).FaceColor = rand(1,3); end;
+legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
+xlabel('Sim. Time'); ylabel('Step Comp. Time [\%]')
+ylim([0,100]); xlim([Ts0D(2),Ts0D(end)]);
+hold on
+yyaxis right
+p2 = plot(Ts0D(2:end),diff(total_Tc),'--r','LineWidth',1.0);
+ylabel('Step Comp. Time [s]')
+ylim([0,1.1*max(diff(total_Tc))])
+set(gca,'ycolor','r')
+FIGNAME = 'profiler';
+save_figure
+end
+
+if 0
+ %% Histograms
+fig = figure;
+histogram(diff(rhs_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(adv_field_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(ghost_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(clos_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(coll_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(poisson_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(Sapj_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(checkfield_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(diag_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+grid on;
+legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
+xlabel('Step Comp. Time [s]'); ylabel('')
+FIGNAME = 'profiler';
+save_figure
+end
\ No newline at end of file
diff --git a/matlab/write_sbash_daint.m b/matlab/write_sbash_daint.m
index 0425cf4c..98c1e359 100644
--- a/matlab/write_sbash_daint.m
+++ b/matlab/write_sbash_daint.m
@@ -48,7 +48,8 @@ fprintf(fid,[...
'module load cray-mpich\n',...
'module load craype-x86-skylake\n',...
'module load cray-fftw\n',...
-'srun ./../../../bin/helaz']);
+'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]);
+%'srun ./../../../bin/helaz']);
fclose(fid);
system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index dd7db865..b4fdf01b 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -35,10 +35,7 @@ fprintf(fid,[...
'#SBATCH --output=out.txt\n',...
'#SBATCH --account=FUA35_TSVVT421\n',...
'#SBATCH --partition=skl_fua_',CLUSTER.PART,'\n',...
-'module load intel\n',...
-'module load intelmpi\n',...
-'module load autoload hdf5/1.10.4--intelmpi--2018--binary\n',...
-'module load fftw\n',...
+'module load autoload hdf5 fftw\n',...
'srun --cpu-bind=cores ./../../../bin/helaz ',num2str(NP_P),' ',num2str(NP_KR)]);
fclose(fid);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 50769e03..4eca5eac 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,25 +1,23 @@
%% Load results
outfile ='';
-if 1
+if 0
%% Load from Marconi
-outfile ='';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_20_J_3_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.8_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.8_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.7_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.7_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_SGGK_CLOS_0_mu_2e-02/out.txt';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_eta_0.6_nu_1e-01/200x100_L_120_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_80_P_10_J_5_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_4e-03/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e+00/200x100_L_120_P_12_J_6_eta_0.6_nu_1e+00_DGGK_CLOS_0_mu_2e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_1e-02/out.txt';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
BASIC.RESDIR = load_marconi(outfile);
end
if 0
%% Load from Daint
- outfile ='';
+ outfile ='/scratch/snx3000/ahoffman/HeLaZ/results/HeLaZ_v2.5_eta_0.6_nu_1e-01/200x100_L_120_P_12_J_6_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_2e-02/out.txt';
BASIC.RESDIR = load_daint(outfile);
end
%%
-% JOBNUM = 0; load_results;
% JOBNUM = 1; load_results;
compile_results
-load_params
%% Retrieving max polynomial degree and sampling info
Npe = numel(Pe); Nje = numel(Je); [JE,PE] = meshgrid(Je,Pe);
@@ -109,7 +107,10 @@ PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz
epsilon_e_pj = zeros(Npe,Nje,Ns5D);
epsilon_i_pj = zeros(Npi,Nji,Ns5D);
-phi_max = zeros(1,Ns2D); % Time evol. of the norm of phi
+phi_maxr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
+phi_avgr_maxz = zeros(1,Ns2D); % Time evol. of the norm of phi
+phi_maxr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
+phi_avgr_avgz = zeros(1,Ns2D); % Time evol. of the norm of phi
Ne_norm = zeros(Npe,Nje,Ns5D); % Time evol. of the norm of Napj
Ni_norm = zeros(Npi,Nji,Ns5D); % .
@@ -117,7 +118,10 @@ Ddr = 1i*KR; Ddz = 1i*KZ; lapl = Ddr.^2 + Ddz.^2;
for it = 1:numel(Ts2D) % Loop over 2D arrays
NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
- phi_max(it) = max(max(squeeze(phi(:,:,it))));
+ phi_maxr_maxz(it) = max( max(squeeze(phi(:,:,it))));
+ phi_avgr_maxz(it) = max(mean(squeeze(phi(:,:,it))));
+ phi_maxr_avgz(it) = mean( max(squeeze(phi(:,:,it))));
+ phi_avgr_avgz(it) = mean(mean(squeeze(phi(:,:,it))));
ExB(it) = max(max(max(abs(phi(3:end,:,it)-phi(1:end-2,:,it))/(2*dr))),max(max(abs(phi(:,3:end,it)-phi(:,1:end-2,it))'/(2*dz))));
GFlux_ri(it) = sum(sum(ni00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_zi(it) = sum(sum(-ni00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
@@ -138,8 +142,9 @@ end
%% Compute growth rate
disp('- growth rate')
% Find max value of transport (end of linear mode)
-[~,itmax] = max(GFlux_ri);
-tstart = 0.1 * Ts2D(itmax); tend = 0.9 * Ts2D(itmax);
+[tmp,tmax] = max(GGAMMA_RI*(2*pi/Nr/Nz)^2);
+[~,itmax] = min(abs(Ts2D-tmax));
+tstart = 0.1 * Ts2D(itmax); tend = 0.5 * Ts2D(itmax);
g_ = zeros(Nkr,Nkz);
for ikr = 1:Nkr
for ikz = 1:Nkz
@@ -158,6 +163,9 @@ if 1
%% Time evolutions and growth rate
fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM),'_',PARAMS];
set(gcf, 'Position', [100, 100, 900, 800])
+subplot(111);
+ suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
+ ', $L=',num2str(L),'$, $N=',num2str(Nr),'$, $(P,J)=(',num2str(PMAXI),',',num2str(JMAXI),')$']);
subplot(421);
for ip = 1:Npe
for ij = 1:Nje
@@ -189,16 +197,20 @@ set(gcf, 'Position', [100, 100, 900, 800])
% plot(Ts5D,PFLUX_RI,'--');
legend(['Gyro. flux';'Part. flux']);
grid on; xlabel('$t c_s/R$'); ylabel('$\Gamma_{r,i}$')
+% ylim([0,2.0]);
subplot(223)
- plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on;
- grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show');
+ plot(kz,g_(1,:),'-','DisplayName','Linear growth rate'); hold on;
+ plot([max(kz)*2/3,max(kz)*2/3],[0,10],'--k', 'DisplayName','2/3 Orszag AA');
+ grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); legend('show');
+% ylim([0,max(g_(1,:))]); xlim([0,max(kz)]);
subplot(224)
- plotname = '$\max_{r,z}(\phi)(t)$';
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
- plot(Ts2D,phi_max,'DisplayName',plotname); hold on;
- grid on; xlabel('$t c_s/R$'); ylabel('$\max_{r,z}(\phi)$'); %legend('show');
-suptitle(['$\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB)]);
+ plot(Ts2D,phi_maxr_maxz,'DisplayName','$\max_{r,z}(\phi)$'); hold on;
+ plot(Ts2D,phi_maxr_avgz,'DisplayName','$\max_{r}\langle\phi\rangle_z$'); hold on;
+ plot(Ts2D,phi_avgr_maxz,'DisplayName','$\max_{z}\langle\phi\rangle_r$'); hold on;
+ plot(Ts2D,phi_avgr_avgz,'DisplayName','$\langle\phi\rangle_{r,z}$'); hold on;
+ grid on; xlabel('$t c_s/R$'); ylabel('$T_e/e$'); legend('show');
save_figure
end
@@ -311,7 +323,7 @@ end
if 0
%% Hermite energy spectra
% tf = Ts2D(end-3);
-fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1000, 300]);
+fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 1800, 600]);
plt = @(x) squeeze(x);
for ij = 1:Nji
subplotnum = 100+Nji*10+ij;
@@ -348,14 +360,14 @@ end
%%
no_AA = (2:floor(2*Nkr/3));
-tKHI = 200;
+tKHI = 100;
[~,itKHI] = min(abs(Ts2D-tKHI));
after_KHI = (itKHI:Ns2D);
if 0
%% Phi frequency space time diagram at kz=0
fig = figure; FIGNAME = ['phi_freq_diag_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
[TY,TX] = meshgrid(Ts2D(after_KHI),kr(no_AA));
- pclr = pcolor(TX,TY,log10(squeeze(abs(PHI(no_AA,1,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar;
+ pclr = pcolor(TX,TY,(squeeze(abs(PHI(no_AA,1,(after_KHI)))))); set(pclr, 'edgecolor','none'); colorbar;
ylabel('$t c_s/R$'), xlabel('$0<k_r<2/3 k_r^{\max}$')
legend('$\log|\tilde\phi(k_z=0)|$')
title('Spectrogram of $\phi$')
@@ -365,7 +377,7 @@ t0 = 0;
[~, it02D] = min(abs(Ts2D-t0));
[~, it05D] = min(abs(Ts5D-t0));
skip_ = 2;
-DELAY = 0.01*skip_;
+DELAY = 0.005*skip_;
FRAMES_2D = it02D:skip_:numel(Ts2D);
FRAMES_5D = it05D:skip_:numel(Ts5D);
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -376,7 +388,7 @@ FIELD = real(ni00); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$n_i$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
-if 1
+if 0
%% Phi real space
GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 1;
FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
@@ -398,12 +410,12 @@ FIELDNAME = '$|\tilde\phi|$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
create_gif
end
if 0
-%% phi @ z = 0
-GIFNAME = ['phi_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
-FIELD =(squeeze(real(phi(:,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
+%% phi averaged on z
+GIFNAME = ['phi_z0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING=1;
+FIELD =(squeeze(mean(real(phi),2))); linestyle = '-.k'; FRAMES = FRAMES_2D;
X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
-FIELDNAME = '$\phi(r=0)$'; XNAME = '$r/\rho_s$';
-create_gif_1D
+FIELDNAME = '$\langle\phi\rangle_{z}$'; XNAME = '$r/\rho_s$';
+create_gif_1D_phi
end
if 0
%% phi @ kz = 0
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
index 44f96b17..24841850 100644
--- a/wk/compute_collision_mat.m
+++ b/wk/compute_collision_mat.m
@@ -29,16 +29,17 @@ kperpmax = sqrt(2) * kmax;
%%
n_ = 1;
for k_ = kperp
+ disp(['----------Computing matrix ',num2str(n_),'/',num2str(Nperp),'----------'])
%% Script to run COSOlver in order to create needed collision matrices
COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
- COSOLVER.pmaxe = 10;
- COSOLVER.jmaxe = 5;
- COSOLVER.pmaxi = 10;
- COSOLVER.jmaxi = 5;
+ COSOLVER.pmaxe = 20;
+ COSOLVER.jmaxe = 10;
+ COSOLVER.pmaxi = 20;
+ COSOLVER.jmaxi = 10;
COSOLVER.kperp = k_;
COSOLVER.neFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2))); % rule of thumb for sum truncation
- COSOLVER.niFLR = max(5,ceil(COSOLVER.kperp^2));
+ COSOLVER.niFLR = min(ceil((2/3*kperpmax)^2),max(5,ceil(COSOLVER.kperp^2)));
COSOLVER.idxT4max = 40;
COSOLVER.neFLRs = 0; % ... only for GK abel
@@ -83,4 +84,5 @@ for k_ = kperp
disp('..done');
cd ../../../HeLaZ/wk
end
+ n_ = n_ + 1;
end
\ No newline at end of file
diff --git a/wk/daint_run.m b/wk/daint_run.m
index 772458f0..e7f53ac3 100644
--- a/wk/daint_run.m
+++ b/wk/daint_run.m
@@ -5,42 +5,53 @@ addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLUSTER PARAMETERS
CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
-CLUSTER.NODES = '1'; % number of nodes
-CLUSTER.NTPN = '36'; % N tasks per node (mpi processes)
-CLUSTER.NTPC = '1'; % N tasks per core (openmp threads)
-CLUSTER.CPUPT = '1'; % CPU per task (number of CPU per mpi proc)
+NP_P = 2; % MPI processes along p
+NP_KR = 18; % MPI processes along kr
CLUSTER.PART = 'normal'; % debug or normal
if(strcmp(CLUSTER.PART,'debug')); CLUSTER.TIME = '00:30:00'; end;
CLUSTER.MEM = '12GB'; % Memory
CLUSTER.JNAME = 'gamma_inf';% Job name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
-ETAB = 0.7; % Magnetic gradient
-NU_HYP = 0.1; % Hyperdiffusivity coefficient
+ETAB = 0.6; % Magnetic gradient
+NU_HYP = 1.0; % Hyperdiffusivity coefficient
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 120; % Size of the squared frequency domain
-P = 10; % Electron and Ion highest Hermite polynomial degree
-J = 05; % Electron and Ion highest Laguerre polynomial degree
+L = 120; % Size of the squared frequency domain
+P = 12; % Electron and Ion highest Hermite polynomial degree
+J = 06; % Electron and Ion highest Laguerre polynomial degree
MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 2500; % Maximal time unit
+TMAX = 1000; % Maximal time unit
DT = 1e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1/2; % Sampling per time unit for 2D arrays
-SPS5D = 1/10; % Sampling per time unit for 5D arrays
-SPSCP = 0; % Sampling per time unit for checkpoints
+SPS0D = 1; % Sampling per time unit for profiler
+SPS2D = 1; % Sampling per time unit for 2D arrays
+SPS5D = 1/40; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
RESTART = 1; % To restart from last checkpoint
-JOB2LOAD= 1;
+JOB2LOAD= 2;
%% OPTIONS
-% SIMID = ['debug']; % Name of the simulation
-SIMID = ['Daint_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
+SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
+% SIMID = 'test_marconi_sugama'; % Name of the simulation
SIMID = sprintf(SIMID,NU);
-CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+PREFIX =[];
+% PREFIX = sprintf('%d_%d_',NP_P, NP_KR);
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
+CO = 1;
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
+NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
+INIT_PHI= 1; % Start simulation with a noisy phi and moments
+%% OUTPUTS
+W_DOUBLE = 1;
+W_GAMMA = 1;
+W_PHI = 1;
+W_NA00 = 1;
+W_NAPJ = 1;
+W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% fixed parameters (for current study)
@@ -60,8 +71,19 @@ NOISE0 = 1.0e-5;
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;
+Nnodes = ceil(Ntot/36);
+Nppn = Ntot/Nnodes;
+CLUSTER.NODES = num2str(Nnodes); % MPI process along p
+CLUSTER.NTPN = num2str(Nppn); % MPI process along kr
+CLUSTER.CPUPT = '1'; % CPU per task
+CLUSTER.NTPC = '1'; % N tasks per core (openmp threads)
%% Run file management scripts
setup
write_sbash_daint
system('rm fort.90 setup_and_run.sh batch_script.sh');
disp('done');
+if(mod(NP_P*NP_KR,36)~= 0)
+ disp('WARNING : unused cores (ntot cores must be a 36 multiple)');
+end
\ No newline at end of file
diff --git a/wk/local_run.m b/wk/local_run.m
index a529cd2f..6812cba4 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -5,33 +5,37 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 1.0; % Collision frequency
-ETAB = 0.5; % Magnetic gradient
+ETAB = 0.6; % Magnetic gradient
ETAN = 1.0; % Density gradient
NU_HYP = 1.0;
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 120; % Size of the squared frequency domain
-PMAXE = 10; % Highest electron Hermite polynomial degree
-JMAXE = 05; % Highest '' Laguerre ''
-PMAXI = 10; % Highest ion Hermite polynomial degree
-JMAXI = 05; % Highest '' Laguerre ''
+L = 80; % Size of the squared frequency domain
+PMAXE = 02; % Highest electron Hermite polynomial degree
+JMAXE = 01; % Highest '' Laguerre ''
+PMAXI = 02; % Highest ion Hermite polynomial degree
+JMAXI = 01; % Highest '' Laguerre ''
+MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0.0/JMAXI^2; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 50; % Maximal time unit
+TMAX = 2500; % Maximal time unit
DT = 2e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
SPS2D = 1/2; % Sampling per time unit for 2D arrays
SPS5D = 1/4; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints/10
-RESTART = 0; % To restart from last checkpoint
-JOB2LOAD= 0;
+RESTART = 1; % To restart from last checkpoint
+JOB2LOAD= 3;
%% OPTIONS AND NAMING
% Collision operator
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 2 ;
+CO = 1;
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
NL_CLOS = -1; % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-SIMID = 'test_SGGK_full_size'; % Name of the simulation
-NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
+SIMID = 'test_stability'; % Name of the simulation
+% SIMID = ['local_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
+% SIMID = sprintf(SIMID,NU);
+NON_LIN = 1; % activate non-linearity (is cancelled if KREQ0 = 1)
%% OUTPUTS
W_DOUBLE = 0;
W_GAMMA = 1;
@@ -47,15 +51,13 @@ KR0KH = 0; A0KH = 0; % Background phi mode to drive Ray-Tay inst.
KREQ0 = 0; % put kr = 0
KPAR = 0.0; % Parellel wave vector component
LAMBDAD = 0.0;
-kmax = N*pi/L;% Highest fourier mode
+kmax = 2/3*N*pi/L;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
% kmaxcut = 2.5;
MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
NOISE0 = 1.0e-5;
TAU = 1.0; % e/i temperature ratio
ETAT = 0.0; % Temperature gradient
-MU_P = 0.0/PMAXI^2; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
-MU_J = 0.0/JMAXI^3; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
INIT_PHI= 1; % Start simulation with a noisy phi and moments
%% Setup and file management
setup
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 450f2d11..346e96e5 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -4,12 +4,13 @@ addpath(genpath('../matlab')) % ... add
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CLUSTER PARAMETERS
-CLUSTER.TIME = '00:10:00'; % allocation time hh:mm:ss
-CLUSTER.PART = 'dbg'; % dbg or prod
-CLUSTER.MEM = '16GB'; % Memory
+CLUSTER.PART = 'prod'; % dbg or prod
+CLUSTER.TIME = '24:00:00'; % allocation time hh:mm:ss
+if(strcmp(CLUSTER.PART,'dbg')); CLUSTER.TIME = '00:30:00'; end;
+CLUSTER.MEM = '64GB'; % Memory
CLUSTER.JNAME = 'gamma_inf';% Job name
-NP_P = 1; % MPI processes along p
-NP_KR = 1; % MPI processes along kr
+NP_P = 2; % MPI processes along p
+NP_KR = 24; % MPI processes along kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
NU = 0.1; % Collision frequency
@@ -17,28 +18,28 @@ ETAB = 0.6; % Magnetic gradient
NU_HYP = 1.0; % Hyperdiffusivity coefficient
%% GRID PARAMETERS
N = 200; % Frequency gridpoints (Nkr = N/2)
-L = 120; % Size of the squared frequency domain
-P = 04; % Electron and Ion highest Hermite polynomial degree
-J = 04; % Electron and Ion highest Laguerre polynomial degree
+L = 80; % Size of the squared frequency domain
+P = 12; % Electron and Ion highest Hermite polynomial degree
+J = 06; % Electron and Ion highest Laguerre polynomial degree
MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 120; % Maximal time unit
-DT = 2e-2; % Time step
+TMAX = 4000; % Maximal time unit
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for profiler
-SPS2D = 1; % Sampling per time unit for 2D arrays
-SPS5D = 1/40; % Sampling per time unit for 5D arrays
+SPS2D = 1/2; % Sampling per time unit for 2D arrays
+SPS5D = 1/100; % Sampling per time unit for 5D arrays
SPSCP = 0; % Sampling per time unit for checkpoints
-RESTART = 0; % To restart from last checkpoint
-JOB2LOAD= 0;
+RESTART = 1; % To restart from last checkpoint
+JOB2LOAD= 4;
%% OPTIONS
-% SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
-SIMID = 'test_marconi_sugama'; % Name of the simulation
+SIMID = ['HeLaZ_v2.5_eta_',num2str(ETAB),'_nu_%0.0e']; % Name of the simulation
SIMID = sprintf(SIMID,NU);
+% SIMID = 'test_marconi_sugama'; % Name of the simulation
PREFIX =[];
% PREFIX = sprintf('%d_%d_',NP_P, NP_KR);
% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
-CO = 2;
+CO = 1;
CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
NL_CLOS = 1; % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
KERN = 0; % Kernel model (0 : GK)
--
GitLab