From 05dc14e817ec2f4bed03882896bbcb1cdbee4eae Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Thu, 1 Apr 2021 13:32:13 +0200
Subject: [PATCH] scripts updates

---
 wk/analysis_2D.m  | 100 +++++++++++++++++++++++-----------------------
 wk/linear_study.m |  53 ++++++++++++++----------
 wk/local_run.m    |  29 ++++++++------
 wk/marconi_run.m  |  25 +++++++-----
 4 files changed, 112 insertions(+), 95 deletions(-)

diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index bd544c99..6d43bb33 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -2,10 +2,10 @@
 outfile ='';
 if 0
     %% Load from Marconi
-    outfile ='';
-    outfile ='';
-    outfile ='';
-    outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_2_12_eta_0.6_nu_1e-01/50x25_L_100_P_20_J_1_eta_0.6_nu_1e-01_DGGK_CLOS_0_mu_3e-01/out.txt';
+outfile ='';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/HeLaZ_v2.4_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.4_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.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';
     BASIC.RESDIR = load_marconi(outfile);
 end
 if 0
@@ -103,6 +103,9 @@ 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
 Ne_norm  = zeros(Npe,Nje,Ns5D);  % Time evol. of the norm of Napj
@@ -124,6 +127,8 @@ 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
@@ -176,7 +181,7 @@ set(gcf, 'Position',  [100, 100, 900, 800])
     end
     grid on; ylabel('$\sum_{k_r,k_z}|N_i^{pj}|$'); xlabel('$t c_s/R$')
     subplot(222)
-        semilogy(Ts0D,GGAMMA_RI*(2*pi/Nr/Nz)^2); hold on;
+        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,'--');
@@ -186,16 +191,10 @@ set(gcf, 'Position',  [100, 100, 900, 800])
         plot(kz,g_(1,:),'-','DisplayName','$\gamma$'); hold on;
         grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma R/c_s$'); %legend('show');
     subplot(224)
-    for ip = 1:Npi
-        for ij = 1:Nji
-            plt      = @(x) squeeze(x(ip,ij,:));
-            plotname = '$\langle\phi\rangle_{r,z}(t)$';
-            clr      = line_colors(min(ip,numel(line_colors(:,1))),:);
-            lstyle   = line_styles(min(ij,numel(line_styles)));
-            plot(Ts2D,phi_max,'DisplayName',plotname,...
-                'Color',clr,'LineStyle',lstyle{1}); hold on;
-        end
-    end
+        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)]);
 save_figure
@@ -247,33 +246,38 @@ end
 
 if 0
 %% Photomaton : real space
-% FIELD = ni00; FNAME = 'ni';
-% FIELD = ne00; FNAME = 'ne';
-FIELD = phi; FNAME = 'phi';
-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));
+% 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((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        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((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        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((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        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((RR),(ZZ),DATA); set(pclr, 'edgecolor','none');pbaspect([1 1 1])
+        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),...
@@ -303,38 +307,34 @@ end
 
 %%
 if 0
-%% Ion moments max mode vs pj
-% tf = Ts2D(end-3); 
-for tf = []
-[~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
-% it2 = it2 + 1;
-fig = figure; FIGNAME = ['kmaxp_Nipj_',sprintf('t=%.2f',Ts2D(it2)),'_',PARAMS];set(gcf, 'Position',  [100, 100, 700, 600])
-
-plt = @(x) squeeze(max(abs(x),[],4));
-% plt = @(x) squeeze(max(fftshift(abs(x),2),[],4));
-
-for ij_ = 1:numel(Ji)
-    subplot(100+numel(Ji)*10+ij_)
-        pclr = imagesc(kr,Pi,plt(Nipj(:,ij_,:,:,it5)));
-        xlabel('$k_r$');
-        if ij_ == 1
-            ylabel('$P$(max o. $k_z$)');
-        else
-            yticks([])
-        end
-        LEGEND = ['$|\hat n_i^{p',num2str(ij_-1),'}|$']; title(LEGEND);
+%% Hermite energy spectra
+% tf = Ts2D(end-3);
+time_array = [1, 100, 400, 1000];
+fig = figure; FIGNAME = ['hermite_spectrum_',PARAMS];set(gcf, 'Position',  [100, 100, 1000, 300]);
+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 |N_i^{p',num2str(Ji(ij)),'}|^2$']; title(TITLE);
 end
 save_figure
 end
-end
 
 
 %%
 t0    = 0;
 [~, it02D] = min(abs(Ts2D-t0));
 [~, it05D] = min(abs(Ts5D-t0));
-skip_ = 1; 
-DELAY = 0.02*skip_;
+skip_ = 2; 
+DELAY = 0.01*skip_;
 FRAMES_2D = it02D:skip_:numel(Ts2D);
 FRAMES_5D = it05D:skip_:numel(Ts5D);
 %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -345,7 +345,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 0
+if 1
 %% Phi real space
 GIFNAME = ['phi',sprintf('_%.2d',JOBNUM),'_',PARAMS];INTERP = 1;
 FIELD = real(phi); X = RR; Y = ZZ; T = Ts2D; FRAMES = FRAMES_2D;
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 06f646f3..9e3b3b6e 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -3,38 +3,49 @@ addpath(genpath('../matlab')) % ... add
 default_plots_options
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
+CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.5;   % Collision frequency
+NU      = 0.1;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.5;
+ETAB    = 0.6;
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
+NU_HYP  = 0.0;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 200;     % Frequency gridpoints (Nkr = N/2)
+N       = 10;     % Frequency gridpoints (Nkr = N/2)
 L       = 120;     % Size of the squared frequency domain
 KREQ0   = 1;      % put kr = 0
 MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARMETERS
 TMAX    = 200;  % Maximal time unit
-DT      = 3e-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   = 1;    % Sampling per time unit for 5D arrays
+SPS5D   = 1/2;    % Sampling per time unit for 5D arrays
 SPSCP   = 0;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 00;
 %% OPTIONS
-SIMID   = 'linear_study_test';  % Name of the simulation
+SIMID   = 'linear_study_SugamaGK';  % Name of the simulation
 NON_LIN = 0 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
-CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+% Collision operator
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
+CO      = 2;
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
+NL_CLOS = 0;   % nonlinear closure model (0: =0 nmax = jmax, 1: nmax = jmax-j, >1 : nmax = NL_CLOS)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_PHI= 0;   % Start simulation with a noisy phi and moments
+%% OUTPUTS
+W_DOUBLE = 0;
+W_GAMMA  = 1;
+W_PHI    = 1;
+W_NA00   = 1;
+W_NAPJ   = 1;
+W_SAPJ   = 0;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % unused
@@ -48,31 +59,31 @@ MU      = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
 %% PARAMETER SCANS
 if 1
 %% Parameter scan over PJ
-% PA = [2, 3, 4, 6, 8, 10];
-% JA = [1, 2, 2, 3, 4,  5];
-% DTA= DT./sqrt(JA);
+PA = [2, 3, 4, 6, 8, 10];
+JA = [1, 2, 2, 3, 4,  5];
+DTA= DT./sqrt(JA);
 mup_ = MU_P;
 muj_ = MU_J;
-PA = [8];
-JA = [4];
+% PA = [4];
+% JA = [1];
 Nparam = numel(PA);
 param_name = 'PJ';
-gamma_Ni00 = zeros(Nparam,N/2+1);
-gamma_Ni21 = zeros(Nparam,N/2+1);
+gamma_Ni00 = zeros(Nparam,floor(N/2)+1);
+gamma_Ni21 = zeros(Nparam,floor(N/2)+1);
 Bohm_transport = zeros(Nparam,1);
-Ni00_ST  = zeros(Nparam,N/2+1,TMAX);
+Ni00_ST  = zeros(Nparam,floor(N/2)+1,SPS2D*TMAX);
 for i = 1:Nparam
     % Change scan parameter
     PMAXE = PA(i); PMAXI = PA(i);
     JMAXE = JA(i); JMAXI = JA(i);
-%     DT = DTA(i);
+    DT = DTA(i);
     MU_P = mup_/PMAXE^2;
     MU_J = muj_/JMAXE^3;
     setup
     % Run linear simulation
-    system(...
-        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk']...
-    )
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz 1 1; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 2 3; cd ../../../wk'])
     % Load and process results
     load_results
     tend   = Ts2D(end); tstart   = 0.4*tend;
@@ -94,7 +105,7 @@ end
 
 if 1
 %% Plot
-SCALE = sqrt(2);
+SCALE = 1;%sqrt(2);
 fig = figure; FIGNAME = 'linear_study';
 plt = @(x) x;
 subplot(211)
diff --git a/wk/local_run.m b/wk/local_run.m
index 9d897ea7..aa822462 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,19 +4,19 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.5;   % Collision frequency
-ETAB    = 0.6;    % Magnetic gradient
+NU      = 0.01;   % Collision frequency
+ETAB    = 0.8;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
-NU_HYP  = 0.1;
+NU_HYP  = 1.0;
 %% GRID PARAMETERS
 N       = 50;     % Frequency gridpoints (Nkr = N/2)
 L       = 100;     % Size of the squared frequency domain
-PMAXE   = 10;     % Highest electron Hermite polynomial degree
-JMAXE   = 1;     % Highest ''       Laguerre ''
-PMAXI   = 10;     % Highest ion      Hermite polynomial degree
-JMAXI   = 1;     % Highest ''       Laguerre ''
+PMAXE   = 4;     % Highest electron Hermite polynomial degree
+JMAXE   = 4;     % Highest ''       Laguerre ''
+PMAXI   = 4;     % Highest ion      Hermite polynomial degree
+JMAXI   = 4;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS
-TMAX    = 100;  % Maximal time unit
+TMAX    = 2000;  % 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
@@ -25,14 +25,17 @@ SPSCP   = 0;    % Sampling per time unit for checkpoints/10
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS AND NAMING
-% SIMID   = ['local_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
-% SIMID   = sprintf(SIMID,NU);
-% SIMID   = 'test_init_phi';  % Name of the simulation
-SIMID   = 'test_parallel_p';  % Name of the simulation
-CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+% Collision operator
+% (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 (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
 KERN    = 0;   % Kernel model (0 : GK)
 INIT_PHI= 1;   % Start simulation with a noisy phi and moments
+% SIMID   = ['local_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
+% SIMID   = sprintf(SIMID,NU);
+% SIMID   = 'test_init_phi';  % Name of the simulation
+SIMID   = ['test_nonlin_NL_CLOS_',num2str(NL_CLOS)];  % Name of the simulation
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 1;
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 575bb761..78b30966 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -4,41 +4,44 @@ addpath(genpath('../matlab')) % ... add
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% CLUSTER PARAMETERS
-CLUSTER.TIME  = '20:00:00'; % allocation time hh:mm:ss
+CLUSTER.TIME  = '12:00:00'; % allocation time hh:mm:ss
 CLUSTER.PART  = 'prod';     % dbg or prod
 CLUSTER.MEM   = '16GB';     % Memory
 CLUSTER.JNAME = 'gamma_inf';% Job name
-NP_P          = 2;          % MPI processes along p  
+NP_P          = 1;          % 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  = 0.1;   % Hyperdiffusivity coefficient
+NU_HYP  = 10.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
+P       = 04;    % Electron and Ion highest Hermite polynomial degree
+J       = 04;    % 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    = 500;  % Maximal time unit
+TMAX    = 2000;  % Maximal time unit
 DT      = 1e-2;  % Time step
 SPS0D   = 1;      % Sampling per time unit for profiler
-SPS2D   = 1/10;   % Sampling per time unit for 2D arrays
+SPS2D   = 1;   % Sampling per time unit for 2D arrays
 SPS5D   = 1/50;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;  % Sampling per time unit for checkpoints
 RESTART = 0;     % To restart from last checkpoint
-JOB2LOAD= 0;
+JOB2LOAD= 1;
 %% OPTIONS
 SIMID   = ['HeLaZ_v2.4_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
 % SIMID   = 'Marconi_parallel_scaling_2D';  % Name of the simulation
 SIMID   = sprintf(SIMID,NU);
 PREFIX  =[];
 % PREFIX  = sprintf('%d_%d_',NP_P, NP_KR);
-CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+% Collision operator
+% (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
+CO      = -3;  
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
+NL_CLOS = 0;   % 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
@@ -69,7 +72,7 @@ ETAN    = 1.0;    % Density gradient
 TAU     = 1.0;    % e/i temperature ratio
 % Compute processes distribution
 Ntot = NP_P * NP_KR;
-Nnodes = ceil(Ntot/24);
+Nnodes = ceil(Ntot/48);
 Nppn   = Ntot/Nnodes; 
 CLUSTER.NODES =  num2str(Nnodes);  % MPI process along p
 CLUSTER.NTPN  =  num2str(Nppn); % MPI process along kr
@@ -79,6 +82,6 @@ setup
 write_sbash_marconi
 system('rm fort.90 setup_and_run.sh batch_script.sh');
 disp('done');
-if(mod(NP_P*NP_KR,24)~= 0)
+if(mod(NP_P*NP_KR,48)~= 0)
     disp('WARNING : unused cores (ntot cores must be a 24 multiple)');
 end
\ No newline at end of file
-- 
GitLab