diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 3e4bb91..a2689bc 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -1,89 +1,103 @@
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_ = [];
PGAMMA_ = [];
PHI_ = [];
Ts0D_ = [];
Ts2D_ = [];
Ts5D_ = [];
Sipj_ = []; Sepj_ = [];
Pe_old = 1e9; Pi_old = Pe_old; Je_old = Pe_old; Ji_old = Pe_old;
while(CONTINUE)
filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
if exist(filename, 'file') == 2
% Load results of simulation #JOBNUM
load_results
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
Pi_new= numel(Pi); Ji_new= numel(Ji);
% If a degree is larger than previous job, put them in a larger array
if (sum([Pe_new, Je_new, Pi_new, Ji_new]>[Pe_old, Je_old, Pi_old, Ji_old]) >= 1)
if W_NAPJ
tmp = Nipj_; sz = size(tmp);
Nipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
Nipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
tmp = Nepj_; sz = size(tmp);
Nepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
Nepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
end
if W_SAPJ
tmp = Sipj_; sz = size(tmp);
Sipj_ = zeros(cat(1,[Pi_new,Ji_new]',sz(3:end)')');
Sipj_(1:Pi_old,1:Ji_old,:,:,:) = tmp;
tmp = Sepj_; sz = size(tmp);
Sepj_ = zeros(cat(1,[Pe_new,Je_new]',sz(3:end)')');
Sepj_(1:Pe_old,1:Je_old,:,:,:) = tmp;
end
end
if W_GAMMA
GGAMMA_ = cat(1,GGAMMA_,GGAMMA_RI);
PGAMMA_ = cat(1,PGAMMA_,PGAMMA_RI);
Ts0D_ = cat(1,Ts0D_,Ts0D);
end
if W_PHI || W_NA00
Ts2D_ = cat(1,Ts2D_,Ts2D);
end
if W_PHI
PHI_ = cat(3,PHI_,PHI);
end
if W_NA00
Ni00_ = cat(3,Ni00_,Ni00);
Ne00_ = cat(3,Ne00_,Ne00);
end
if W_NAPJ || W_SAPJ
Ts5D_ = cat(1,Ts5D_,Ts5D);
end
if W_NAPJ
Nipj_ = cat(5,Nipj_,Nipj);
Nepj_ = cat(5,Nepj_,Nepj);
end
if W_SAPJ
Sipj_ = cat(5,Sipj_,Sipj);
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)
CONTINUE = 0;
disp(['found ',num2str(JOBFOUND),' results']);
end
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_;
Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_;
clear Nipj_ Nepj_ Ni00_ Ne00_ PHI_ Ts2D_ Ts5D_ GGAMMA_ PGAMMA_ Ts0D_
Sipj = Sipj_; Sepj = Sepj_;
clear Sipj_ Sepj_
JOBNUM = LASTJOB
filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
\ No newline at end of file
diff --git a/matlab/create_gif_1D_phi.m b/matlab/create_gif_1D_phi.m
new file mode 100644
index 0000000..f68761c
--- /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 37d363b..57b8470 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -1,61 +1,95 @@
%% load profiling
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
CPUTIME = double(h5readatt(filename,'/data/input','cpu_time'));
DT_SIM = h5readatt(filename,'/data/input','dt');
rhs_Tc = h5read(filename,'/profiler/Tc_rhs');
adv_field_Tc = h5read(filename,'/profiler/Tc_adv_field');
ghost_Tc = h5read(filename,'/profiler/Tc_ghost');
clos_Tc = h5read(filename,'/profiler/Tc_clos');
coll_Tc = h5read(filename,'/profiler/Tc_coll');
poisson_Tc = h5read(filename,'/profiler/Tc_poisson');
Sapj_Tc = h5read(filename,'/profiler/Tc_Sapj');
checkfield_Tc= h5read(filename,'/profiler/Tc_checkfield');
diag_Tc = h5read(filename,'/profiler/Tc_diag');
step_Tc = h5read(filename,'/profiler/Tc_step');
Ts0D = h5read(filename,'/profiler/time');
missing_Tc = step_Tc - rhs_Tc - adv_field_Tc - ghost_Tc -clos_Tc ...
-coll_Tc -poisson_Tc -Sapj_Tc -checkfield_Tc -diag_Tc;
total_Tc = step_Tc;
TIME_PER_FCT = [diff(rhs_Tc); diff(adv_field_Tc); diff(ghost_Tc);...
diff(clos_Tc); diff(coll_Tc); diff(poisson_Tc); diff(Sapj_Tc); ...
diff(checkfield_Tc); diff(diag_Tc); diff(missing_Tc)];
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))
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 0425cf4..98c1e35 100644
--- a/matlab/write_sbash_daint.m
+++ b/matlab/write_sbash_daint.m
@@ -1,56 +1,57 @@
% Write the input script "fort.90" with desired parameters
INPUT = 'setup_and_run.sh';
fid = fopen(INPUT,'wt');
fprintf(fid,[...
'#!/bin/bash\n',...
'mkdir -p $SCRATCH/HeLaZ/wk\n',...
...
'cd $SCRATCH/HeLaZ/wk/\n',...
...
'mkdir -p ', BASIC.RESDIR,'\n',...
'cd ',BASIC.RESDIR,'\n',...
'cp $HOME/HeLaZ/wk/fort.90 .\n',...
'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
...
'jid=$(sbatch batch_script.sh)\n',...
'echo $jid\n',...
'echo to check output log :\n',...
'echo tail -f $SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
fclose(fid);
system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
% Write the sbatch script
INPUT = 'batch_script.sh';
fid = fopen(INPUT,'wt');
fprintf(fid,[...
'#!/bin/bash\n',...
'#SBATCH --job-name="',CLUSTER.JNAME,'"\n',...
'#SBATCH --time=', CLUSTER.TIME,'\n',...
'#SBATCH --nodes=', CLUSTER.NODES,'\n',...
'#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
'#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
'#SBATCH --ntasks-per-core=', CLUSTER.NTPC,'\n',...
'#SBATCH --mem=', CLUSTER.MEM,'\n',...
'#SBATCH --error=err.txt\n',...
'#SBATCH --output=out.txt\n',...
'#SBATCH --account="s882"\n',...
'#SBATCH --constraint=mc\n',...
'#SBATCH --hint=nomultithread\n',...
'#SBATCH --partition=',CLUSTER.PART,'\n',...
'export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n',...
...% '#SBATCH --job-name=',PARAMS,'\n\n',...
'module purge\n',...
'module load PrgEnv-intel\n',...
'module load cray-hdf5-parallel\n',...
'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,'/.']);
system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@ela.cscs.ch:HeLaZ/wk');
\ No newline at end of file
diff --git a/matlab/write_sbash_marconi.m b/matlab/write_sbash_marconi.m
index dd7db86..b4fdf01 100644
--- a/matlab/write_sbash_marconi.m
+++ b/matlab/write_sbash_marconi.m
@@ -1,47 +1,44 @@
% Write the input script "fort.90" with desired parameters
INPUT = 'setup_and_run.sh';
fid = fopen(INPUT,'wt');
fprintf(fid,[...
'#!/bin/bash\n',...
'mkdir -p $CINECA_SCRATCH/HeLaZ/wk\n',...
...
'cd $CINECA_SCRATCH/HeLaZ/wk/\n',...
...
'mkdir -p ', BASIC.RESDIR,'\n',...
'cd ',BASIC.RESDIR,'\n',...
'cp $HOME/HeLaZ/wk/fort.90 .\n',...
'cp $HOME/HeLaZ/wk/batch_script.sh .\n',...
...
'sbatch batch_script.sh\n',...
'echo tail -f $CINECA_SCRATCH/HeLaZ/results/',BASIC.SIMID,'/',BASIC.PARAMS,'/out.txt']);
fclose(fid);
system(['cp setup_and_run.sh ',BASIC.RESDIR,'/.']);
% Write the sbatch script
INPUT = 'batch_script.sh';
fid = fopen(INPUT,'wt');
fprintf(fid,[...
'#!/bin/bash\n',...
'#SBATCH --job-name=',CLUSTER.JNAME,'\n',...
'#SBATCH --time=', CLUSTER.TIME,'\n',...
'#SBATCH --nodes=', CLUSTER.NODES,'\n',...
'#SBATCH --cpus-per-task=', CLUSTER.CPUPT,'\n',...
'#SBATCH --ntasks-per-node=', CLUSTER.NTPN,'\n',...
'#SBATCH --mem=', CLUSTER.MEM,'\n',...
'#SBATCH --error=err.txt\n',...
'#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);
system(['cp batch_script.sh ',BASIC.RESDIR,'/.']);
system('scp {fort.90,setup_and_run.sh,batch_script.sh} ahoffman@login.marconi.cineca.it:/marconi/home/userexternal/ahoffman/HeLaZ/wk');
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 50769e0..4eca5ea 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,454 +1,466 @@
%% 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);
Npi = numel(Pi); Nji = numel(Ji); [JI,PI] = meshgrid(Ji,Pi);
Ns5D = numel(Ts5D);
Ns2D = numel(Ts2D);
% renaming and reshaping quantity of interest
Ts5D = Ts5D';
Ts2D = Ts2D';
%% Build grids
Nkr = numel(kr); Nkz = numel(kz);
[KZ,KR] = meshgrid(kz,kr);
Lkr = max(kr)-min(kr); Lkz = max(kz)-min(kz);
dkr = Lkr/(Nkr-1); dkz = Lkz/(Nkz-1);
KPERP2 = KZ.^2+KR.^2;
[~,ikr0] = min(abs(kr)); [~,ikz0] = min(abs(kz));
Lk = max(Lkr,Lkz);
dr = 2*pi/Lk; dz = 2*pi/Lk;
Nr = max(Nkr,Nkz); Nz = Nr;
r = dr*(-Nr/2:(Nr/2-1)); Lr = max(r)-min(r);
z = dz*(-Nz/2:(Nz/2-1)); Lz = max(z)-min(z);
[ZZ,RR] = meshgrid(z,r);
%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Analysis :')
disp('- iFFT')
% IFFT (Lower case = real space, upper case = frequency space)
ne00 = zeros(Nr,Nz,Ns2D); % Gyrocenter density
ni00 = zeros(Nr,Nz,Ns2D);
np_i = zeros(Nr,Nz,Ns5D); % Ion particle density
si00 = zeros(Nr,Nz,Ns5D);
phi = zeros(Nr,Nz,Ns2D);
drphi = zeros(Nr,Nz,Ns2D);
dr2phi = zeros(Nr,Nz,Ns2D);
dzphi = zeros(Nr,Nz,Ns2D);
for it = 1:numel(Ts2D)
NE_ = Ne00(:,:,it); NI_ = Ni00(:,:,it); PH_ = PHI(:,:,it);
ne00(:,:,it) = real(fftshift(ifft2((NE_),Nr,Nz)));
ni00(:,:,it) = real(fftshift(ifft2((NI_),Nr,Nz)));
phi (:,:,it) = real(fftshift(ifft2((PH_),Nr,Nz)));
drphi(:,:,it) = real(fftshift(ifft2(1i*KR.*(PH_),Nr,Nz)));
dr2phi(:,:,it)= real(fftshift(ifft2(-KR.^2.*(PH_),Nr,Nz)));
dzphi(:,:,it) = real(fftshift(ifft2(1i*KZ.*(PH_),Nr,Nz)));
end
% Building a version of phi only 5D sampling times
PHI_Ts5D = zeros(Nkr,Nkz,Ns5D);
err = 0;
for it = 1:numel(Ts5D) % Loop over 5D arrays
[shift, it2D] = min(abs(Ts2D-Ts5D(it)));
if shift > err; err = shift; end;
PHI_Ts5D(:,:,it) = PHI(:,:,it2D);
end
if err > 0; disp('WARNING Ts2D and Ts5D are shifted'); end;
Np_i = zeros(Nkr,Nkz,Ns5D); % Ion particle density in Fourier space
for it = 1:numel(Ts5D)
[~, it2D] = min(abs(Ts2D-Ts5D(it)));
Np_i(:,:,it) = 0;
for ij = 1:Nji
Kn = (KPERP2/2.).^(ij-1) .* exp(-KPERP2/2)/(factorial(ij-1));
Np_i(:,:,it) = Np_i(:,:,it) + Kn.*squeeze(Nipj(1,ij,:,:,it));
end
np_i(:,:,it) = real(fftshift(ifft2(squeeze(Np_i(:,:,it)),Nr,Nz)));
end
% Post processing
disp('- post processing')
E_pot = zeros(1,Ns2D); % Potential energy n^2
E_kin = zeros(1,Ns2D); % Kinetic energy grad(phi)^2
ExB = zeros(1,Ns2D); % ExB drift intensity \propto |\grad \phi|
% gyrocenter and particle flux from real space
GFlux_ri = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni drphi>
GFlux_zi = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ni dzphi>
GFlux_re = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne drphi>
GFlux_ze = zeros(1,Ns2D); % Gyrocenter flux Gamma = <ne dzphi>
PFlux_ri = zeros(1,Ns5D); % Particle flux
% gyrocenter and particle flux from fourier coefficients
GFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Ni00.*conj(PHI),1),2)))*(2*pi/Nr/Nz)^2;
PFLUX_RI = real(squeeze(sum(sum(-1i*KZ.*Np_i.*conj(PHI_Ts5D),1),2)))*(2*pi/Nr/Nz)^2;
% Hermite energy spectrum
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); % .
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;
GFlux_re(it) = sum(sum(ne00(:,:,it).*dzphi(:,:,it)))*dr*dz/Lr/Lz;
GFlux_ze(it) = sum(sum(-ne00(:,:,it).*drphi(:,:,it)))*dr*dz/Lr/Lz;
end
for it = 1:numel(Ts5D) % Loop over 5D arrays
[~, it2D] = min(abs(Ts2D-Ts5D(it)));
Ne_norm(:,:,it)= sum(sum(abs(Nepj(:,:,:,:,it)),3),4)/Nkr/Nkz;
Ni_norm(:,:,it)= sum(sum(abs(Nipj(:,:,:,:,it)),3),4)/Nkr/Nkz;
epsilon_e_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nepj(:,:,:,:,it)).^2,3),4);
epsilon_i_pj(:,:,it) = sqrt(pi)/2*sum(sum(abs(Nipj(:,:,:,:,it)).^2,3),4);
% Particle flux
PFlux_ri(it) = sum(sum(np_i(:,:,it).*dzphi(:,:,it2D)))*dr*dz/Lr/Lz;
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
g_(ikr,ikz) = LinearFit_s(Ts2D,squeeze(abs(Ni00(ikr,ikz,:))),tstart,tend);
end
end
[gmax,ikzmax] = max(g_(1,:));
kzmax = abs(kz(ikzmax));
Bohm_transport = ETAB/ETAN*gmax/kzmax^2;
%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
default_plots_options
disp('Plots')
FMT = '.fig';
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
plt = @(x) squeeze(x(ip,ij,:));
plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$'];
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname,...
'Color',clr,'LineStyle',lstyle{1}); hold on;
end
end
grid on; ylabel('$\sum_{k_r,k_z}|N_e^{pj}|$');
subplot(423)
for ip = 1:Npi
for ij = 1:Nji
plt = @(x) squeeze(x(ip,ij,:));
plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$'];
clr = line_colors(min(ip,numel(line_colors(:,1))),:);
lstyle = line_styles(min(ij,numel(line_styles)));
semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname,...
'Color',clr,'LineStyle',lstyle{1}); hold on;
end
end
grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
subplot(222)
plot(Ts0D,GGAMMA_RI*(2*pi/Nr/Nz)^2); hold on;
% plot(Ts2D,GFLUX_RI)
plot(Ts0D,PGAMMA_RI*(2*pi/Nr/Nz)^2);
% 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
if 0
%% Particle fluxes
SCALING = Nkr*dkr * Nkz*dkz;
fig = figure; FIGNAME = ['gamma',sprintf('_%.2d',JOBNUM),'_',PARAMS];
set(gcf, 'Position', [100, 100, 800, 300])
semilogy(Ts2D,GFLUX_RI, 'color', line_colors(2,:)); hold on
plot(Ts5D,PFLUX_RI,'.', 'color', line_colors(2,:)); hold on
plot(Ts2D,SCALING*GFlux_ri, 'color', line_colors(1,:)); hold on
plot(Ts5D,SCALING*PFlux_ri,'.', 'color', line_colors(1,:)); hold on
xlabel('$tc_{s0}/\rho_s$'); ylabel('$\Gamma_r$'); grid on
title(['$\eta=',num2str(ETAB),'\quad',...
'\nu_{',CONAME,'}=',num2str(NU),'$', ' $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'])
legend('Gyro Flux','Particle flux', 'iFFT GFlux', 'iFFT PFlux')%'$\eta\gamma_{max}/k_{max}^2$')
save_figure
end
if 0
%% Space time diagramm (fig 11 Ivanov 2020)
fig = figure; FIGNAME = ['space_time_drphi','_',PARAMS];set(gcf, 'Position', [100, 100, 1200, 600])
subplot(311)
semilogy(Ts2D,GFLUX_RI,'-'); hold on
plot(Ts5D,PFLUX_RI,'.'); hold on
% plot(Ts2D,Bohm_transport*ones(size(Ts2D)),'--'); hold on
ylabel('$\Gamma_r$'); grid on
title(['$\eta=',num2str(ETAB),'\quad',...
'\nu_{',CONAME,'}=',num2str(NU),'$'])
legend(['$P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$'],'Particle flux')%'$\eta\gamma_{max}/k_{max}^2$')
% set(gca,'xticklabel',[])
subplot(312)
yyaxis left
plot(Ts2D,squeeze(max(max((phi)))))
ylabel('$\max \phi$')
yyaxis right
plot(Ts2D,squeeze(mean(max(dr2phi))))
ylabel('$s\sim\langle\partial_r^2\phi\rangle_z$'); grid on
% set(gca,'xticklabel',[])
subplot(313)
[TY,TX] = meshgrid(r,Ts2D);
pclr = pcolor(TX,TY,squeeze(mean(drphi(:,:,:),2))'); set(pclr, 'edgecolor','none'); %colorbar;
xlabel('$t c_s/R$'), ylabel('$r/\rho_s$')
legend('$\langle\partial_r \phi\rangle_z$')
save_figure
end
if 0
%% Photomaton : real space
% FIELD = ni00; FNAME = 'ni'; XX = RR; YY = ZZ;
FIELD = phi; FNAME = 'phi'; XX = RR; YY = ZZ;
% FIELD = fftshift(abs(Ni00),2); FNAME = 'Fni'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
% FIELD = fftshift(abs(PHI),2); FNAME = 'Fphi'; XX = fftshift(KR,2); YY = fftshift(KZ,2);
tf = 100; [~,it1] = min(abs(Ts2D-tf));
tf = 118; [~,it2] = min(abs(Ts2D-tf));
tf = 140; [~,it3] = min(abs(Ts2D-tf));
tf = 300; [~,it4] = min(abs(Ts2D-tf));
fig = figure; FIGNAME = [FNAME,'_snaps','_',PARAMS]; set(gcf, 'Position', [100, 100, 1500, 400])
plt = @(x) x;%./max(max(x));
subplot(141)
DATA = plt(FIELD(:,:,it1));
pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
colormap gray
xlabel('$r/\rho_s$'); ylabel('$z/\rho_s$');set(gca,'ytick',[]);
title(sprintf('$t c_s/R=%.0f$',Ts2D(it1)));
subplot(142)
DATA = plt(FIELD(:,:,it2));
pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
colormap gray
xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
title(sprintf('$t c_s/R=%.0f$',Ts2D(it2)));
subplot(143)
DATA = plt(FIELD(:,:,it3));
pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
colormap gray
xlabel('$r/\rho_s$');ylabel('$z/\rho_s$');set(gca,'ytick',[]);
title(sprintf('$t c_s/R=%.0f$',Ts2D(it3)));
subplot(144)
DATA = plt(FIELD(:,:,it4));
pclr = pcolor((XX),(YY),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
colormap gray
xlabel('$r/\rho_s$');ylabel('$z/\rho_s$'); set(gca,'ytick',[]);
title(sprintf('$t c_s/R=%.0f$',Ts2D(it4)));
% suptitle(['$\',FNAME,'$, $\nu_{',CONAME,'}=$', num2str(NU), ', $\eta_B=$',num2str(ETAB),...
% ', $P=',num2str(PMAXI),'$, $J=',num2str(JMAXI),'$']);
save_figure
end
%%
if 0
%% Show frame in kspace
tf = 0; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
fig = figure; FIGNAME = ['krkz_',sprintf('t=%.0f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position', [100, 100, 700, 600])
subplot(221); plt = @(x) fftshift((abs(x)),2);
pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(PHI(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('$t c_s/R=%.0f$',Ts2D(it2))); legend('$|\hat\phi|$');
subplot(222); plt = @(x) fftshift(abs(x),2);
pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Ni00(:,:,it2))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{00}|$');
subplot(223); plt = @(x) fftshift((abs(x)),2); FIELD = squeeze(Nipj(1,2,:,:,:));
pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(FIELD(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$k_r$'); ylabel('$k_z$'); legend('$|\hat n_i^{pj=01}|$');
subplot(224); plt = @(x) fftshift((abs(x)),2);
pclr = pcolor(fftshift(KR,2),fftshift(KZ,2),plt(Si00(:,:,it5))); set(pclr, 'edgecolor','none'); colorbar;
xlabel('$k_r$'); ylabel('$k_z$');legend('$\hat S_i^{00}$');
save_figure
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;
subplot(subplotnum)
for it5 = 1:2:Ns5D
alpha = it5*1.0/Ns5D;
loglog(Pi(1:2:end),plt(epsilon_i_pj(1:2:end,ij,it5)),...
'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
end
grid on;
xlabel('$p$');
TITLE = ['$\sum_{kr,kz} |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE);
end
save_figure
end
%%
if 0
%% Laguerre energy spectra
% tf = Ts2D(end-3);
fig = figure; FIGNAME = ['laguerre_spectrum_',PARAMS];set(gcf, 'Position', [100, 100, 500, 400]);
plt = @(x) squeeze(x);
for it5 = 1:2:Ns5D
alpha = it5*1.0/Ns5D;
loglog(Ji,plt(max(epsilon_i_pj(:,:,it5),[],1)),...
'color',(1-alpha)*[0.8500, 0.3250, 0.0980]+alpha*[0, 0.4470, 0.7410],...
'DisplayName',['t=',num2str(Ts5D(it5))]); hold on;
end
grid on;
xlabel('$j$');
TITLE = ['$\max_p\sum_{kr,kz} |N_i^{pj}|^2$']; title(TITLE);
save_figure
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$')
end
%%
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% Density ion
GIFNAME = ['ni',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
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;
FIELDNAME = '$\phi$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
%% radial particle transport
GIFNAME = ['gamma_r',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 1;
FIELD = real(ni00.*dzphi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
FIELDNAME = '$\Gamma_r$'; XNAME = '$r/\rho_s$'; YNAME = '$z/\rho_s$';
create_gif
end
if 0
%% Phi fourier
GIFNAME = ['FFT_phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 0;
FIELD = ifftshift((abs(PHI)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D; FRAMES = FRAMES_2D;
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
GIFNAME = ['phi_kz0',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; SCALING = 0;
FIELD =squeeze(log10(abs(PHI(no_AA,1,:)))); linestyle = '-.'; FRAMES = FRAMES_2D;
X = kr(no_AA); T = Ts2D; YMIN = -30; YMAX = 6; XMIN = min(kr); XMAX = max(kr);
FIELDNAME = '$|\tilde\phi(k_z=0)|$'; XNAME = '$k_r\rho_s$';
create_gif_1D
end
if 0
%% Density ion frequency
GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D;
FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
create_gif
end
if 0
%% Density electron frequency
GIFNAME = ['Ne00',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0; FRAMES = FRAMES_2D;
FIELD =ifftshift((abs(Ne00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
FIELDNAME = '$N_e^{00}$'; XNAME = '$k_r\rho_s$'; YNAME = '$k_z\rho_s$';
create_gif
end
if 0
%% kr vs P Si
GIFNAME = ['Sip0_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
plt = @(x) squeeze(max((abs(x)),[],4));
FIELD =plt(Sipj(:,1,:,:,:)); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$';
create_gif_imagesc
end
if 0
%% maxkz, kr vs p, for all Nipj over time
GIFNAME = ['Nipj_kr',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
plt = @(x) squeeze(max((abs(x)),[],4));
FIELD = plt(Nipj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$, ${k_z}^{max}$';
create_gif_5D
end
if 0
%% maxkr, kz vs p, for all Nipj over time
GIFNAME = ['Nipj_kz',sprintf('_%.2d',JOBNUM),'_',PARAMS]; INTERP = 0;
plt = @(x) fftshift(squeeze(max((abs(x)),[],3)),3);
FIELD = plt(Nipj); X = sort(kz'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, ${k_r}^{max}$';
create_gif_5D
end
%%
\ No newline at end of file
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
index 44f96b1..2484185 100644
--- a/wk/compute_collision_mat.m
+++ b/wk/compute_collision_mat.m
@@ -1,86 +1,88 @@
addpath(genpath('../matlab')) % ... add
%% Grid configuration
N = 200; % Frequency gridpoints (Nkr = N/2)
L = 120; % Size of the squared frequency domain
dk = 2*pi/L;
kmax = N/2*dk;
kr = dk*(0:N/2);
kz = dk*(0:N/2);
[KZ, KR]= meshgrid(kz,kr);
KPERP = sqrt(KR.^2 + KZ.^2);
kperp = reshape(KPERP,[1,numel(kr)^2]);
kperp = uniquetol(kperp,1e-14);
Nperp = numel(kperp);
if 0
%% plot the kperp distribution
figure
plot(kperp)
end
% %% Check if the differences btw kperp is larger than naming precision
% dkperp = diff(kperp);
% warning = sum(dkperp<0.0002);
% if warning > 0
% disp('Warning : dkperp < 0.0002');
% end
%%
%% We compute only on a kperp grid with dk space from 0 to kperpmax
kperp = unique([0:dk:(sqrt(2)*kmax),sqrt(2)*kmax]);
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
COSOLVER.npeFLR = 0; % ... only for GK abel
COSOLVER.niFLRs = 0; % ... only for GK abel
COSOLVER.npiFLR = 0; % ... only for GK abel
COSOLVER.gk = 1;
COSOLVER.co = 3;
% ! 0 = nothing
% ! 1 = coulomb
% ! 2 = pitch-angle (with constant coll.freq.)
% ! 3 = sugama
% ! 4 = pitch-angle with veloty dependent coll. freq.
% ! 5 = improved sugama
% ! 6 = Hirschman-Sigmar Clarke
% ! 7 = Abel (only for like species)
% ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
% ! 9 = GK coulomb polarization term
write_fort90_COSOlver
k_string = sprintf('%0.4f',k_);
mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
'_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
'_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
'_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
'_JE_12',...
'_NFLR_',num2str(COSOLVER.neFLR),...
'_kperp_',k_string,'.h5'];
if (exist(mat_file_name,'file') > 0)
disp(['Matrix available for kperp = ',k_string]);
else
cd ../../Documents/MoliSolver/COSOlver/
disp(['Matrix not found for kperp = ',k_string]);
disp([num2str(n_),'/',Nperp])
disp('computing...');
CMD = 'mpirun -np 6 bin/CO 2 2 2 > out.txt';
disp(CMD);
system(CMD);
system(CMD);
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 772458f..e7f53ac 100644
--- a/wk/daint_run.m
+++ b/wk/daint_run.m
@@ -1,67 +1,89 @@
%clear all;
addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 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)
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;
NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
PMAXE = P; % Highest electron Hermite polynomial degree
JMAXE = J; % Highest '' Laguerre ''
PMAXI = P; % Highest ion Hermite polynomial degree
JMAXI = J; % Highest '' Laguerre ''
kmax = N*pi/L;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
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 a529cd2..6812cba 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -1,62 +1,64 @@
addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
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;
W_PHI = 1;
W_NA00 = 1;
W_NAPJ = 1;
W_SAPJ = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% unused
KERN = 0; % Kernel model (0 : GK)
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
system('rm fort.90');
\ No newline at end of file
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 450f2d1..346e96e 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -1,86 +1,87 @@
%clear all;
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
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)
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)
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;
NON_LIN = 1 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
PMAXE = P; % Highest electron Hermite polynomial degree
JMAXE = J; % Highest '' Laguerre ''
PMAXI = P; % Highest ion Hermite polynomial degree
JMAXI = J; % Highest '' Laguerre ''
kmax = N*pi/L;% Highest fourier mode
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
MU = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
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/48);
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
%% Run file management scripts
setup
write_sbash_marconi
system('rm fort.90 setup_and_run.sh batch_script.sh');
disp('done');
if(mod(NP_P*NP_KR,48)~= 0)
disp('WARNING : unused cores (ntot cores must be a 48 multiple)');
end
\ No newline at end of file