From dd62fee629cfff072ac40f4d6e5704afb0fff019 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 15 Mar 2021 16:10:44 +0100
Subject: [PATCH] Added an option to start with only noise in phi for smoother
 starts

---
 matlab/setup.m        |  4 ++--
 matlab/write_fort90.m |  3 ++-
 wk/daint_run.m        | 12 ++++++------
 wk/linear_study.m     | 43 ++++++++++++++++++++++---------------------
 wk/local_run.m        | 23 ++++++++++++-----------
 wk/marconi_run.m      | 21 +++++++++++----------
 6 files changed, 55 insertions(+), 51 deletions(-)

diff --git a/matlab/setup.m b/matlab/setup.m
index 6460fe51..6f42c2b2 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -34,10 +34,10 @@ MODEL.eta_n   = ETAN;        % source term kappa for HW
 MODEL.eta_T   = ETAT;        % Temperature
 MODEL.eta_B   = ETAB;        % Magnetic
 MODEL.lambdaD = LAMBDAD;
-if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
+% if A0KH ~= 0; SIMID = [SIMID,'_Nz_',num2str(L/2/pi*KR0KH),'_A_',num2str(A0KH)]; end;
 % Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
-INITIAL.only_Na00         = '.false.';
+if INIT_PHI; INITIAL.init_noisy_phi = '.true.'; else; INITIAL.init_noisy_phi = '.false.';end;
 INITIAL.initback_moments  = 0.0e-5;
 INITIAL.initnoise_moments = NOISE0;
 INITIAL.iseed             = 42;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index bf34207b..1e2f98da 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -57,7 +57,8 @@ fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&INITIAL_CON\n');
-fprintf(fid,['  only_Na00         =', INITIAL.only_Na00,'\n']);
+fprintf(fid,['  INIT_NOISY_PHI    =', INITIAL.init_noisy_phi,'\n']);
+% fprintf(fid,['  only_Na00         =', '.false.','\n']);
 fprintf(fid,['  initback_moments  =', num2str(INITIAL.initback_moments),'\n']);
 fprintf(fid,['  initnoise_moments =', num2str(INITIAL.initnoise_moments),'\n']);
 fprintf(fid,['  iseed             =', num2str(INITIAL.iseed),'\n']);
diff --git a/wk/daint_run.m b/wk/daint_run.m
index 73efd98c..772458f0 100644
--- a/wk/daint_run.m
+++ b/wk/daint_run.m
@@ -6,7 +6,7 @@ addpath(genpath('../matlab')) % ... add
 %% CLUSTER PARAMETERS
 CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 CLUSTER.NODES = '1';        % number of nodes
-CLUSTER.NTPN  = '32';       % N tasks per node (mpi processes)
+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)
 CLUSTER.PART  = 'normal';   % debug or normal
@@ -15,8 +15,8 @@ CLUSTER.MEM   = '12GB';     % Memory
 CLUSTER.JNAME = 'gamma_inf';% Job name
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
-ETAB    = 0.6;   % Magnetic gradient
+NU      = 0.1;   % Collision frequency
+ETAB    = 0.7;   % Magnetic gradient
 NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
 %% GRID PARAMETERS
 N       = 200;   % Frequency gridpoints (Nkr = N/2)
@@ -26,14 +26,14 @@ J       = 05;    % 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    = 1000;  % Maximal time unit
+TMAX    = 2500;  % 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
-RESTART = 0;     % To restart from last checkpoint
-JOB2LOAD= 0;
+RESTART = 1;     % To restart from last checkpoint
+JOB2LOAD= 1;
 %% OPTIONS
 % SIMID   = ['debug'];  % Name of the simulation
 SIMID   = ['Daint_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 51d8c938..06f646f3 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -5,23 +5,23 @@ default_plots_options
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
+NU      = 0.5;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.6;
+ETAB    = 0.5;
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
 LAMBDAD = 0.0;
 NOISE0  = 1.0e-5;
 %% GRID PARAMETERS
-N       = 150;     % Frequency gridpoints (Nkr = N/2)
-L       = 70;     % Size of the squared frequency domain
+N       = 200;     % 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    = 100;  % Maximal time unit
-DT      = 1e-2;   % Time step
+TMAX    = 200;  % Maximal time unit
+DT      = 3e-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
@@ -29,12 +29,12 @@ SPSCP   = 0;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 00;
 %% OPTIONS
-SIMID   = 'linear_study_test_mu_kin';  % Name of the simulation
+SIMID   = 'linear_study_test';  % 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)
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 KERN    = 0;   % Kernel model (0 : GK)
-
+INIT_PHI= 0;   % Start simulation with a noisy phi and moments
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % unused
@@ -48,13 +48,13 @@ 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 = [4];
-% JA = [2];
+PA = [8];
+JA = [4];
 Nparam = numel(PA);
 param_name = 'PJ';
 gamma_Ni00 = zeros(Nparam,N/2+1);
@@ -65,13 +65,13 @@ 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 4 ./../../../bin/helaz; cd ../../../wk']...
+        ['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz 1 6; cd ../../../wk']...
     )
     % Load and process results
     load_results
@@ -86,40 +86,41 @@ for i = 1:Nparam
     end
     gamma_Ni00(i,:) = real(gamma_Ni00(i,:) .* (gamma_Ni00(i,:)>=0.0));
     gamma_Ni21(i,:) = real(gamma_Ni21(i,:) .* (gamma_Ni21(i,:)>=0.0));
-    kzmax = abs(kr(ikzmax));
-    Bohm_transport(i) = ETAB/ETAN*gmax/kzmax^2;
+%     kzmax = abs(kr(ikzmax));
+%     Bohm_transport(i) = ETAB/ETAN*gmax/kzmax^2;
     % Clean output
     system(['rm -r ',BASIC.RESDIR])
 end
 
 if 1
 %% Plot
+SCALE = sqrt(2);
 fig = figure; FIGNAME = 'linear_study';
 plt = @(x) x;
 subplot(211)
     for i = 1:Nparam
         clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
         linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
-        plot(plt(kr),plt(gamma_Ni00(i,:)),...
+        plot(plt(SCALE*kr),plt(gamma_Ni00(i,:)),...
             'Color',clr,...
             'LineStyle',linestyle{1},...
             'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
         hold on;
     end
-    grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma(N_i^{00})\rho_s/c_s$'); xlim([0.0,max(kr)]);
+    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{00})\rho_s/c_s$'); xlim([0.0,max(kr)]);
     title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
     legend('show')
 subplot(212)
     for i = 1:Nparam
         clr       = line_colors(mod(i-1,numel(line_colors(:,1)))+1,:);
         linestyle = line_styles(floor((i-1)/numel(line_colors(:,1)))+1);
-        plot(plt(kr),plt(gamma_Ni21(i,:)),...
+        plot(plt(SCALE*kr),plt(gamma_Ni21(i,:)),...
             'Color',clr,...
             'LineStyle',linestyle{1},...
             'DisplayName',['$P=$',num2str(PA(i)),', $J=$',num2str(JA(i))]);
         hold on;
     end
-    grid on; xlabel('$k_z\rho_s$'); ylabel('$\gamma(N_i^{21})\rho_s/c_s$'); xlim([0.0,max(kr)]);
+    grid on; xlabel('$k_z\rho_s^{R}$'); ylabel('$\gamma(N_i^{21})\rho_s/c_s$'); xlim([0.0,max(kr)]);
     title(['$\eta_B=',num2str(ETAB),'$, $\nu_{',CONAME,'}=',num2str(NU),'$, ', CLOSNAME])
     legend('show')
 saveas(fig,[SIMDIR,'gamma_Ni_vs_',param_name,'_',PARAMS,'.fig']);
diff --git a/wk/local_run.m b/wk/local_run.m
index 856d8c64..142d0b92 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,34 +4,35 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
-ETAB    = 0.5;    % Magnetic gradient
+NU      = 0.5;   % Collision frequency
+ETAB    = 0.6;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 NU_HYP  = 0.1;
 %% GRID PARAMETERS
-N       = 200;     % Frequency gridpoints (Nkr = N/2)
-L       = 125;     % Size of the squared frequency domain
-PMAXE   = 2;     % Highest electron Hermite polynomial degree
+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   = 2;     % Highest ion      Hermite polynomial degree
+PMAXI   = 10;     % Highest ion      Hermite polynomial degree
 JMAXI   = 1;     % Highest ''       Laguerre ''
-
 %% TIME PARAMETERS
-TMAX    = 50;  % Maximal time unit
-DT      = 1e-2;   % Time step
+TMAX    = 10;  % Maximal time unit
+DT      = 2e-2;   % Time step
 SPS0D   = 1/DT;    % 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
+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;
 %% OPTIONS
 % 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)
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 KERN    = 0;   % Kernel model (0 : GK)
+INIT_PHI= 1;   % Start simulation with a noisy phi and moments
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% unused
@@ -43,7 +44,7 @@ NON_LIN = 1 *(1-KREQ0);   % activate non-linearity (is cancelled if KREQ0 = 1)
 kmax    = N*pi/L;% Highest fourier mode
 HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
 % kmaxcut = 2.5;
-MU      = NU_HYP/(kmax)^4 % Hyperdiffusivity coefficient
+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
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index d53f5473..0b8a24c3 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -4,7 +4,7 @@ addpath(genpath('../matlab')) % ... add
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% CLUSTER PARAMETERS
-CLUSTER.TIME  = '12:00:00'; % allocation time hh:mm:ss
+CLUSTER.TIME  = '24:00:00'; % allocation time hh:mm:ss
 CLUSTER.NODES = '1';        % MPI process
 CLUSTER.CPUPT = '1';        % CPU per task
 CLUSTER.NTPN  = '24';       % N tasks per node
@@ -13,14 +13,14 @@ CLUSTER.MEM   = '16GB';     % Memory
 CLUSTER.JNAME = 'gamma_inf';% Job name
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 1.0;   % Collision frequency
+NU      = 0.1;   % Collision frequency
 ETAB    = 0.6;   % Magnetic gradient
-NU_HYP  = 0.2;   % Hyperdiffusivity coefficient
+NU_HYP  = 0.1;   % Hyperdiffusivity coefficient
 %% GRID PARAMETERS
-N       = 200;     % Frequency gridpoints (Nkr = N/2)
-L       = 70;     % Size of the squared frequency domain
-P       = 10;       % Electron and Ion highest Hermite polynomial degree
-J       = 5;       % Electron and Ion highest Laguerre polynomial degree
+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
 MU_P    = 0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
 MU_J    = 0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% TIME PARAMETERS
@@ -28,17 +28,18 @@ 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/4;  % Sampling per time unit for 5D arrays
+SPS5D   = 1/10;  % Sampling per time unit for 5D arrays
 SPSCP   = 0;  % Sampling per time unit for checkpoints
 RESTART = 0;     % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = ['Marconi_DGGK_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
-% SIMID   = 'Marconi_restart';  % Name of the simulation
+SIMID   = ['Marconi_new_eta_',num2str(ETAB),'_nu_%0.0e'];  % Name of the simulation
+% SIMID   = 'Marconi_test';  % Name of the simulation
 SIMID   = sprintf(SIMID,NU);
 CO      = -3;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
 KERN    = 0;   % Kernel model (0 : GK)
+INIT_PHI= 0;   % Start simulation with a noisy phi and moments
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% fixed parameters (for current study)
-- 
GitLab