From 9c460ecf6b6e14f09695c912a492b07ebf6981e6 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 15 Jan 2021 14:32:20 +0100
Subject: [PATCH] update scripts

---
 Makefile                 |  6 +++-
 matlab/compile_results.m | 25 ++++++++++++++++-
 matlab/load_params.m     | 11 ++++++++
 matlab/setup.m           | 19 +++++++++----
 matlab/write_fort90.m    |  2 +-
 wk/analysis_2D.m         | 59 ++++++++++++++++++----------------------
 wk/linear_study.m        | 58 ++++++++++++---------------------------
 wk/local_run.m           | 22 ++++++++-------
 wk/marconi_run.m         |  8 +++---
 9 files changed, 114 insertions(+), 96 deletions(-)

diff --git a/Makefile b/Makefile
index 456be934..a1e15cad 100644
--- a/Makefile
+++ b/Makefile
@@ -42,7 +42,8 @@ cleanbin:
 $(OBJDIR)/diagnose.o : src/srcinfo.h
 
 FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
-$(OBJDIR)/coeff_mod.o $(OBJDIR)/collision_mod.o $(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
+$(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/collision_mod.o \
+$(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
 $(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
 $(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/main.o $(OBJDIR)/memory.o \
 $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o \
@@ -68,6 +69,9 @@ $(OBJDIR)/utility_mod.o
  $(OBJDIR)/coeff_mod.o : src/coeff_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o
 		$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/coeff_mod.F90 -o $@
 
+ $(OBJDIR)/closure_mod.o : src/closure_mod.F90 $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/fields_mod.o
+	 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/closure_mod.F90 -o $@
+
  $(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o
 	 	$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
 
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index ee40fafe..9b1f9d83 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -6,12 +6,33 @@ PHI_     = [];
 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
+        sz = size(Nepj); Pe_new= sz(1); Je_new= sz(2);
+        sz = size(Nipj); Pi_new= sz(1); Ji_new= sz(2);
+        % 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)
+            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;
+            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
+        
+        
         Nipj_ = cat(5,Nipj_,Nipj);
         Nepj_ = cat(5,Nepj_,Nepj);
         Ni00_ = cat(3,Ni00_,Ni00);
@@ -28,6 +49,8 @@ while(CONTINUE)
         CONTINUE = 0;
         disp(['found ',num2str(JOBNUM),' results']);
     end
+    Pe_old = Pe_new; Je_old = Je_new;
+    Pi_old = Pi_new; Ji_old = Ji_new;
 end
 Nipj = Nipj_; Nepj = Nepj_; Ts5D = Ts5D_;
 Ni00 = Ni00_; Ne00 = Ne00_; PHI = PHI_; Ts2D = Ts2D_;
diff --git a/matlab/load_params.m b/matlab/load_params.m
index bec47edf..09cee32f 100644
--- a/matlab/load_params.m
+++ b/matlab/load_params.m
@@ -22,3 +22,14 @@ if NON_LIN == 'y'
 else
     NON_LIN = 0;
 end
+
+degngrad   = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
+    '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI),...
+    '_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_nu_%0.0e_',...
+    CONAME,'_mu_%0.0e'];
+degngrad   = sprintf(degngrad,[NU,MU]);
+if ~NON_LIN; degngrad = ['lin_',degngrad]; end
+resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
+gridname   = ['L_',num2str(L),'_'];
+PARAMS = [resolution,gridname,degngrad];
+BASIC.RESDIR = [SIMDIR,PARAMS,'/'];
diff --git a/matlab/setup.m b/matlab/setup.m
index dd598a39..f85f571d 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -13,7 +13,7 @@ GRID.kpar  = KPAR;
 
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-if 0;      MODEL.DK      = '.true.'; else; MODEL.DK      = '.false.';end;
+MODEL.CLOS   = CLOS;
 if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
 MODEL.mu      = MU;
 MODEL.nu      = NU; % hyper diffusive coefficient nu for HW
@@ -61,10 +61,19 @@ elseif(CO == -1); CONAME = 'FC';
 elseif(CO == -2); CONAME = 'DG';
 elseif(CO == -3); CONAME = 'DGGK';
 end
-degngrad   = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
-    '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI),...
-    '_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_nu_%0.0e_',...
-    CONAME,'_mu_%0.0e'];
+if    (CLOS == 0); CLOSNAME = 'Trunc.';
+elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
+elseif(CLOS == 2); CLOSNAME = 'Clos. 2';
+end
+if (PMAXE == PMAXI) && (JMAXE == JMAXI)
+    degngrad   = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE)];
+else
+    degngrad   = ['Pe_',num2str(PMAXE),'_Je_',num2str(JMAXE),...
+        '_Pi_',num2str(PMAXI),'_Ji_',num2str(JMAXI)];
+end
+degngrad = [degngrad,'_eta_',num2str(ETAB/ETAN),'_nu_%0.0e_',...
+        CONAME,'_CLOS_',num2str(CLOS),'_mu_%0.0e'];
+    
 degngrad   = sprintf(degngrad,[NU,MU]);
 if ~NON_LIN; degngrad = ['lin_',degngrad]; end
 resolution = [num2str(GRID.Nr),'x',num2str(GRID.Nz/2),'_'];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 23c331ca..d17d1ec0 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -38,7 +38,7 @@ fprintf(fid,'/\n');
 fprintf(fid,'&MODEL_PAR\n');
 fprintf(fid,'  ! Collisionality\n');
 fprintf(fid,['  CO      = ', num2str(MODEL.CO),'\n']);
-fprintf(fid,['  DK      = ', num2str(MODEL.DK),'\n']);
+fprintf(fid,['  CLOS   = ', num2str(MODEL.CLOS),'\n']);
 fprintf(fid,['  NON_LIN = ', MODEL.NON_LIN,'\n']);
 fprintf(fid,['  mu      = ', num2str(MODEL.mu),'\n']);
 fprintf(fid,['  nu      = ', num2str(MODEL.nu),'\n']);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 65c87e85..26c8b526 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -6,8 +6,7 @@ if 0
     outfile ='';
     outfile ='';
     outfile ='';
-    outfile ='';
-    outfile ='';
+    outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/Marconi_DGGK/200x100_L_100_Pe_8_Je_4_Pi_8_Ji_4_nB_0.66_nN_1_nu_1e-01_DGGK_mu_1e-03/out.txt';
     BASIC.RESDIR = load_marconi(outfile);
 end
 %%
@@ -150,7 +149,10 @@ set(gcf, 'Position',  [100, 100, 900, 800])
         for ij = 1:Nje
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$N_e^{',num2str(Pe(ip)),num2str(Je(ij)),'}$'];
-            semilogy(Ts5D,plt(Ne_norm),'DisplayName',plotname); hold on;
+            clr      = line_colors(ip,:);
+            lstyle   = line_styles(ij);
+            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}|$');
@@ -159,7 +161,10 @@ set(gcf, 'Position',  [100, 100, 900, 800])
         for ij = 1:Nji
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$N_i^{',num2str(Pi(ip)),num2str(Ji(ij)),'}$'];
-            semilogy(Ts5D,plt(Ni_norm),'DisplayName',plotname); hold on;
+             clr      = line_colors(ip,:);
+            lstyle   = line_styles(ij);
+            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}|$');
@@ -174,7 +179,10 @@ set(gcf, 'Position',  [100, 100, 900, 800])
         for ij = 1:Nji
             plt      = @(x) squeeze(x(ip,ij,:));
             plotname = ['$S_i^{',num2str(ip-1),num2str(ij-1),'}$'];
-            semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname); hold on;
+            clr      = line_colors(ip,:);
+            lstyle   = line_styles(ij);
+            semilogy(Ts5D,plt(Si_norm),'DisplayName',plotname,...
+                'Color',clr,'LineStyle',lstyle{1}); hold on;
         end
     end
     grid on; xlabel('$t c_s/R$'); ylabel('$\sum_{k_r,k_z}|S_i^{pj}|$'); %legend('show');
@@ -187,10 +195,10 @@ if 0
 % FIELD = ni00; FNAME = 'ni';
 % FIELD = ne00; FNAME = 'ne';
 FIELD = phi; FNAME = 'phi';
-tf = 60;  [~,it1] = min(abs(Ts2D-tf));
-tf = 120;  [~,it2] = min(abs(Ts2D-tf)); 
-tf = 250; [~,it3] = min(abs(Ts2D-tf));
-tf = 400; [~,it4] = min(abs(Ts2D-tf));
+tf = 200;  [~,it1] = min(abs(Ts2D-tf));
+tf = 600;  [~,it2] = min(abs(Ts2D-tf)); 
+tf =1000; [~,it3] = min(abs(Ts2D-tf));
+tf =2000; [~,it4] = min(abs(Ts2D-tf));
 fig = figure; FIGNAME = [FNAME,'_snaps']; set(gcf, 'Position',  [100, 100, 1500, 400])
 plt = @(x) x;%./max(max(x));
     subplot(141)
@@ -276,7 +284,7 @@ fig = figure; FIGNAME = 'space_time_drphi';set(gcf, 'Position',  [100, 100, 1200
     subplot(211)
     yyaxis left
     plot(Ts2D,GFlux_ri); hold on
-    plot(Ts5D,PFlux_ri,'o');
+%     plot(Ts5D,PFlux_ri,'o');
     ylabel('$\Gamma_r$'); grid on
     ylim([0,1.1*max(GFlux_ri)]);
     yyaxis right
@@ -295,34 +303,19 @@ end
 
 %%
 if 0
-%% Mode time evolution
-[~,ikr ] = min(abs(kr-dkr));
-[~,ik00] = min(abs(kz));
-[~,idk]  = min(abs(kz-dkz));
-[~,ik50] = min(abs(kz-0.1*max(kz)));
-[~,ik75] = min(abs(kz-0.2*max(kz)));
-[~,ik10] = min(abs(kz-0.3*max(kz)));
-plt = @(x) abs(squeeze(x));
-fig = figure; FIGNAME = ['mode_time_evolution',sprintf('_%.2d',JOBNUM)];
-        semilogy(Ts2D,plt(Ni00(ikr,ik00,:)),'DisplayName', ...
-            ['$k_z = $',num2str(kz(ik00))]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr,idk,:)),'DisplayName', ...
-            ['$k_z = $',num2str(kz(idk))]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr,ik50,:)),'DisplayName', ...
-            ['$k_z = $',num2str(kz(ik50))]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr,ik75,:)),'DisplayName', ...
-            ['$k_z = $',num2str(kz(ik75))]); hold on
-        semilogy(Ts2D,plt(Ni00(ikr,ik10,:)),'DisplayName', ...
-            ['$k_z = $',num2str(kz(ik10))]); hold on
-        xlabel('$t$'); ylabel('$\hat n_i^{00}$'); legend('show');
-title(sprintf('$k_r=$ %1.1f',kr(ikr)))
-save_figure
+%% Space time diagram of moments amplitude
+pj_grid_i = 1:((PMAXI+1)*(JMAXI+1));
+[TY,TX] = meshgrid(Ts5D,pj_grid_i);
+Ni_norm_tmp = squeeze(reshape(Ni_norm,[(PMAXI+1)*(JMAXI+1),1,numel(Ts5D)]));
+fig = figure; FIGNAME = 'space_time_Ni_norm';
+% pclr = image(pj_grid_i,flip(Ts5D),Ni_norm_tmp'); hold on;
+pclr = pcolor(TX,TY,log(Ni_norm_tmp)); hold on; set(pclr, 'edgecolor','none'); %colorbar;
 end
 
 %%
 if 0
 %% Show frame in kspace
-tf = 300; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
+tf = 1000; [~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['krkz_frame',sprintf('t=%.0f',Ts2D(it2))];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;
diff --git a/wk/linear_study.m b/wk/linear_study.m
index f290e719..92ace15e 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -5,7 +5,7 @@ default_plots_options
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e+0;   % Collision frequency
+NU      = 0e-01;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 0.5;
 ETAN    = 1.0;    % Density gradient
@@ -14,12 +14,12 @@ MU      = 1e-3;   % Hyper diffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 100;     % Frequency gridpoints (Nkr = N/2)
-L       = 100;     % Size of the squared frequency domain
+N       = 50;     % Frequency gridpoints (Nkr = N/2)
+L       = 10;     % Size of the squared frequency domain
 KREQ0   = 1;      % put kr = 0
 %% TIME PARMETERS
 TMAX    = 100;  % Maximal time unit
-DT      = 5e-2;   % Time step
+DT      = 1e-2;   % Time step
 SPS0D   = 0.5;      % Sampling per time unit for 2D arrays
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
@@ -29,7 +29,10 @@ JOB2LOAD= 00;
 %% OPTIONS
 SIMID   = 'linear_study';  % Name of the simulation
 NON_LIN = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
-CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+CLOS    = 2;   % Closure model (0 : =0 truncation, 1 : n+j = min(nmax,n+j), 2: odd/even adapted)
+KERN    = 0;   % Kernel model (0 : GK)
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % unused
@@ -40,10 +43,11 @@ KPAR    = 0.0;    % Parellel wave vector component
 %% PARAMETER SCANS
 if 1
 %% Parameter scan over PJ
-% PA = [2, 6, 8, 12];
-% JA = [1, 3, 4,  6];
-PA = [4];
-JA = [2];
+PA = [2, 3, 4, 6, 12, 20];
+JA = [1, 2, 2, 3,  6, 10];
+DTA= DT./sqrt(JA);
+% PA = [4];
+% JA = [2];
 Nparam = numel(PA);
 param_name = 'PJ';
 gamma_Ni = zeros(Nparam,N/2+1);
@@ -52,10 +56,11 @@ for i = 1:Nparam
     % Change scan parameter
     PMAXE = PA(i); PMAXI = PA(i);
     JMAXE = JA(i); JMAXI = JA(i);
+    DT = DTA(i);
     setup
     % Run linear simulation
     system(...
-        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz; cd ../../../wk']...
+        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 4 ./../../../bin/helaz; cd ../../../wk']...
     )
     % Load and process results
     load_results
@@ -66,36 +71,7 @@ for i = 1:Nparam
     end
     gamma_Ni(i,:) = real(gamma_Ni(i,:) .* (gamma_Ni(i,:)>=0.0));
     % Clean output
-    %system(['rm -r ',BASIC.RESDIR])
-end
-% if 0
-% %% Plot
-% fig = figure; FIGNAME = 'space_time_Ni00';
-% i = 1;
-%
-% [YY, XX] = meshgrid(Ts2D,kz);
-% plt = @(x) squeeze(log(abs(x)))
-% pclr = pcolor(XX,YY,plt(Ni00_ST(i,:,:)));set(pclr, 'edgecolor','none');
-% grid on; xlabel('$k_z$'); ylabel('$t$');
-% title(['$\eta_B=',num2str(ETAB),'$,$P,J=',num2str(PA(i)),',',JA(i),'$'])
-% legend('show')
-% saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
-% saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
-% end
-if 0
-%% Plot
-fig = figure; FIGNAME = 'space_time_Ni00';
-i = 3;
-plt = @(x) squeeze((abs(x)));
-for ikr=1:4:numel(kr)/2+1
-    NAME = ['$k_z=',num2str(kr(ikr)),'$'];
-    semilogy(1:TMAX,plt(Ni00_ST(i,ikr,1:TMAX)), 'DisplayName', NAME); hold on;
-end
-grid on; xlabel('$t c_s/\rho_s$'); ylabel('$|Ni^{00}|$');
-title(['$\eta_B=',num2str(ETAB),'$; $P,J=',num2str(PA(i)),',',num2str(JA(i)),'$, $\nu=',num2str(NU),'$'])
-legend('show');
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
-saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
+    system(['rm -r ',BASIC.RESDIR])
 end
 if 1
 %% Plot
@@ -111,7 +87,7 @@ for i = 1:Nparam
     hold on;
 end
 grid on; xlabel('$k_z$'); ylabel('$\gamma(N_i^{00})$'); xlim([0.0,max(kr)]);
-title(['$\eta_B=',num2str(ETAB),'$, $\nu=',num2str(NU),'$'])
+title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
 legend('show')
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.png']);
diff --git a/wk/local_run.m b/wk/local_run.m
index 8a8244a8..61af9e38 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,7 +4,7 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1e+1;   % Collision frequency
+NU      = 1e-1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
 ETAB    = 0.5;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
@@ -14,22 +14,24 @@ NU_HYP  = 0.1;
 %% GRID PARAMETERS
 N       = 128;     % Frequency gridpoints (Nkr = N/2)
 L       = 66;     % Size of the squared frequency domain
-PMAXE   = 3;     % Highest electron Hermite polynomial degree
-JMAXE   = 2;     % Highest ''       Laguerre ''
-PMAXI   = 3;     % Highest ion      Hermite polynomial degree
-JMAXI   = 2;     % Highest ''       Laguerre ''
+PMAXE   = 4;     % Highest electron Hermite polynomial degree
+JMAXE   = 3;     % Highest ''       Laguerre ''
+PMAXI   = 4;     % Highest ion      Hermite polynomial degree
+JMAXI   = 3;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
-TMAX    = 100;  % Maximal time unit
+TMAX    = 500;  % Maximal time unit
 DT      = 3e-2;   % Time step
 SPS0D   = 1/DT;    % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/5;    % Sampling per time unit for 5D arrays
 SPSCP   = 1/10;    % Sampling per time unit for checkpoints
-RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 0;
+RESTART = 1;      % To restart from last checkpoint
+JOB2LOAD= 3;
 %% OPTIONS
-SIMID   = 'test_cleaning';  % Name of the simulation
-CO      = -1;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+% SIMID   = 'test_pp2=p_closure';  % Name of the simulation
+SIMID   = 'test_truncations';  % Name of the simulation
+CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+CLOS   = 0;   % Truncation method (0 : =0 closure, 1 : n+j = min(nmax,n+j))
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index f817c8d2..bc527e45 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -19,8 +19,8 @@ NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
 %% GRID PARAMETERS
 N       = 200;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
-P       = 8;       % Electron and Ion highest Hermite polynomial degree
-J       = 4;       % Electron and Ion highest Laguerre polynomial degree
+P       = 6;       % Electron and Ion highest Hermite polynomial degree
+J       = 3;       % Electron and Ion highest Laguerre polynomial degree
 %% TIME PARAMETERS
 TMAX    = 150;  % Maximal time unit
 DT      = 1e-2;   % Time step
@@ -28,12 +28,12 @@ SPS0D   = 10;    % Sampling per time unit for profiler
 SPS2D   = 1;      % Sampling per time unit for 2D arrays
 SPS5D   = 1/10;    % Sampling per time unit for 5D arrays
 SPSCP   = 1/10;    % Sampling per time unit for checkpoints
-RESTART = 0;      % To restart from last checkpoint
+RESTART = 1;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
 SIMID   = 'Marconi_DGGK';  % Name of the simulation
 CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
-
+CLOS   = 0;   % Truncation method (0 : =0 closure, 1 : n+j = min(nmax,n+j))
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-- 
GitLab