From ab81f3e047fa39e1af8b858a74ff1202aafc773a Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 4 Aug 2021 09:54:14 +0200
Subject: [PATCH] scripts update

---
 matlab/setup.m      | 19 +++++++------
 wk/HD_study.m       |  7 +++++
 wk/analysis_3D.m    | 12 +++++---
 wk/linear_study.m   | 39 +++++++++++++-------------
 wk/local_run.m      |  8 +++---
 wk/marconi_run.m    | 16 +++++------
 wk/plot_cosol_mat.m | 68 +++++++++++++++++++++++++++++++--------------
 7 files changed, 104 insertions(+), 65 deletions(-)

diff --git a/matlab/setup.m b/matlab/setup.m
index 0f6ded1f..821068bd 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -54,19 +54,22 @@ if (abs(CO) == 2) %Sugama operator
 elseif (abs(CO) == 3) %pitch angle operator
     INITIAL.mat_file = ['''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'''];
 elseif (CO == 4) % Full Coulomb GK
-    disp('Warning, FCGK not implemented yet')
+    INITIAL.mat_file = ['''../../../iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5'''];
 elseif (CO == -1) % DGDK
     disp('Warning, DGDK not debugged')
 end
 
 % Naming and creating input file
-if    (CO == -3); CONAME = 'PADK';
-elseif(CO == -2); CONAME = 'SGDK';
-elseif(CO == -1); CONAME = 'DGDK';
-elseif(CO ==  0); CONAME = 'LBDK';
-elseif(CO ==  1); CONAME = 'DGGK';
-elseif(CO ==  2); CONAME = 'SGGK';
-elseif(CO ==  3); CONAME = 'PAGK';
+switch abs(CO)
+    case 0; CONAME = 'LB';
+    case 1; CONAME = 'DG';
+    case 2; CONAME = 'SG';
+    case 3; CONAME = 'PA';
+    case 4; CONAME = 'FC';
+    otherwise; CONAME ='UK';
+end
+if (CO <= 0); CONAME = [CONAME,'DK'];
+else;         CONAME = [CONAME,'GK'];
 end
 if    (CLOS == 0); CLOSNAME = 'Trunc.';
 elseif(CLOS == 1); CLOSNAME = 'Clos. 1';
diff --git a/wk/HD_study.m b/wk/HD_study.m
index dbf3f4e2..7ec342c5 100644
--- a/wk/HD_study.m
+++ b/wk/HD_study.m
@@ -34,6 +34,7 @@ rp_=[...
     5e-2, 1e-3;...
     1e-2, 1e-3;...
     1e-2, 1e-4;...
+    1e-2, 5e-4;...
     ];
 figure; set(gcf, 'Position',  [100, 100, 900, 400])
 grid on; xlim([5e-4,5e0]); ylim([5e-5,5e+0]);
@@ -72,6 +73,12 @@ mu_ = [3e-3 1e-3 1e-4];
 nu_ = [1e-2 1e-2 1e-2];
 plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=2,1');
 
+% HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_CLOS_0_mu_3e-02/
+mu_ = [3e-3 5e-4];
+nu_ = [1e-2 1e-2];
+plot(nu_,mu_,'x--','DisplayName','N=150, L=100, P,J=4,2');
+
+
 % HD_study/150x75_L_100_P_2_J_1_eta_0.6_nu_5e-02_DGGK_CLOS_0_mu_3e-02/
 mu_ = [3e-3 1e-3];
 nu_ = [5e-2 5e-2];
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index 5300342d..083fd6e3 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -7,8 +7,10 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='HD_study/300x150_L_100_P_2_J_1_eta_0.6_nu_5e-01_DGGK_mu_2e-03';
-% outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
@@ -19,7 +21,9 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/test_3D_marconi/100x50_L_60_P_2_J_1_eta_Inf_nu_1e-01_DGGK_mu_0e+00/out.txt';
+outfile ='';
+outfile ='';
+outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/test_3D_marconi/100x50x20_L_60_q0_1_P_2_J_1_eta_0.6_nu_1e-01_DGGK_mu_2e-02/out.txt';
 BASIC.RESDIR      = ['../',outfile(46:end-8),'/'];
 BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
 end
@@ -74,7 +78,7 @@ t0    =00; iz = 1; ix = 1; iy = 1;
 skip_ = 1; DELAY = 1e-2*skip_;
 [~, it03D] = min(abs(Ts3D-t0)); FRAMES_3D = it03D:skip_:numel(Ts3D);
 [~, it05D] = min(abs(Ts5D-t0)); FRAMES_5D = it05D:skip_:numel(Ts5D);
-INTERP = 0; T = Ts3D; FRAMES = FRAMES_3D;
+INTERP = 1; T = Ts3D; FRAMES = FRAMES_3D;
 % Field to plot
 % FIELD = dens_e; NAME = 'ne';   FIELDNAME = 'n_e';
 % FIELD = dens_i; NAME = 'ni';   FIELDNAME = 'n_i';
diff --git a/wk/linear_study.m b/wk/linear_study.m
index 223568a4..81d50007 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -1,6 +1,6 @@
-for NU = [0.01]
+for NU = [0.1]
 for ETAN = [2.0]
-for CO = [3]
+for CO = [1]
 %clear all;
 addpath(genpath('../matlab')) % ... add
 default_plots_options
@@ -33,12 +33,10 @@ SPSCP   = 0;    % Sampling per time unit for checkpoints
 RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 00;
 %% OPTIONS
-% SIMID   = 'v2.7_lin_analysis';  % Name of the simulation
-SIMID   = 'kobayashi_2015_fig2';  % Name of the simulation
-% SIMID   = 'v2.6_lin_analysis';  % Name of the simulation
+SIMID   = 'v3.1_lin_analysis';  % Name of the simulation
 NON_LIN = 0 *(1-KXEQ0);   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % Collision operator
-% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Full Couloumb ; +/- for GK/DK)
+% (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle, 4 : Full Couloumb ; +/- for GK/DK)
 % CO      = 2;
 INIT_ZF = 0; ZF_AMP = 0.0;
 CLOS    = 0;   % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
@@ -63,15 +61,18 @@ JOBNUM  = 00;
 KPAR    = 0.0;    % Parellel wave vector component
 HD_CO   = 0.5;    % Hyper diffusivity cutoff ratio
 kmax    = N*pi/L;% Highest fourier mode
-MU      = NU_HYP/(HD_CO*kmax)^4 % Hyperdiffusivity coefficient
-
+MU      = NU_HYP/(HD_CO*kmax)^4; % Hyperdiffusivity coefficient
+Nz      = 1;      % number of perpendicular planes (parallel grid)
+q0      = 1.0;    % safety factor
+shear   = 0.0;    % magnetic shear
+eps     = 0.0;    % inverse aspect ratio
 %% PARAMETER SCANS
 if 1
 %% Parameter scan over PJ
 % PA = [2 4 6 10];
 % JA = [1 2 3  5];
-PA = [2];
-JA = [1];
+PA = [6];
+JA = [3];
 DTA= DT*ones(size(JA));%./sqrt(JA);
 % DTA= DT;
 mup_ = MU_P;
@@ -87,23 +88,21 @@ for i = 1:Nparam
     PMAXE = PA(i); PMAXI = PA(i);
     JMAXE = JA(i); JMAXI = JA(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 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'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 1 ./../../../bin/helaz_3 1 1; cd ../../../wk'])
+    system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 1 6; cd ../../../wk'])
+%     system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ./../../../bin/helaz_3 2 3; cd ../../../wk'])
 %     Load and process results
     %%
     filename = ['../results/',SIMID,'/',PARAMS,'/outputs_00.h5'];
     load_results
-    tend   = Ts2D(end); tstart   = 0.4*tend;
-    [~,itstart] = min(abs(Ts2D-tstart));
-    [~,itend]   = min(abs(Ts2D-tend));
+    tend   = Ts3D(end); tstart   = 0.4*tend;
+    [~,itstart] = min(abs(Ts3D-tstart));
+    [~,itend]   = min(abs(Ts3D-tend));
     for ikx = 1:N/2+1
-        gamma_Ni00(i,ikx) = (LinearFit_s(Ts2D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend))))));
-        Ni00_ST(i,ikx,1:numel(Ts2D)) = squeeze((Ni00(ikx,1,:)));
+        gamma_Ni00(i,ikx) = (LinearFit_s(Ts3D(itstart:itend)',(squeeze(abs(Ni00(ikx,1,itstart:itend))))));
+        Ni00_ST(i,ikx,1:numel(Ts3D)) = squeeze((Ni00(ikx,1,:)));
     end
     tend   = Ts5D(end); tstart   = 0.4*tend;
     [~,itstart] = min(abs(Ts5D-tstart));
diff --git a/wk/local_run.m b/wk/local_run.m
index c588c465..f6f67373 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -6,10 +6,10 @@ CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %% PHYSICAL PARAMETERS
 NU      = 0.1;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
-NU_HYP  = 0.1;
+NU_HYP  = 1.0;
 %% GRID AND GEOMETRY PARAMETERS
-N       = 100;     % Frequency gridpoints (Nkx = N/2)
-L       = 60;     % Size of the squared frequency domain
+N       = 50;     % Frequency gridpoints (Nkx = N/2)
+L       = 30;     % Size of the squared frequency domain
 Nz      = 10;      % number of perpendicular planes (parallel grid)
 q0      = 1.0;    % safety factor
 shear   = 0.0;    % magnetic shear
@@ -17,7 +17,7 @@ eps     = 0.0;    % inverse aspect ratio
 P       = 2;
 J       = 1;
 %% TIME PARAMETERS
-TMAX    = 200;  % Maximal time unit
+TMAX    = 500;  % 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
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 236dc473..f93e63cf 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -18,13 +18,13 @@ NP_P          = 2;          % MPI processes along p
 NP_KX         = 24;         % MPI processes along kx
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.05;   % Collision frequency
+NU      = 0.1;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
-NU_HYP  = 10.0;
+NU_HYP  = 1.0;
 %% GRID PARAMETERS
-N       = 300;     % Frequency gridpoints (Nkx = N/2)
-L       = 100;     % Size of the squared frequency domain
-Nz      = 1;      % number of perpendicular planes (parallel grid)
+N       = 100;     % Frequency gridpoints (Nkx = N/2)
+L       = 60;     % Size of the squared frequency domain
+Nz      = 20;      % number of perpendicular planes (parallel grid)
 q0      = 1.0;    % q factor ()
 shear   = 0.0;    % magnetic shear
 eps     = 0.0;    % inverse aspect ratio
@@ -43,11 +43,11 @@ JOB2LOAD= 0;
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO      = 2;
+CO      = 1;
 CLOS    = 0;   % Closure model (0: =0 truncation)
 NL_CLOS = 0;   % nonlinear closure model (-2: nmax = jmax, -1: nmax = jmax-j, >=0 : nmax = NL_CLOS)
-% SIMID   = 'test_3D_marconi';  % Name of the simulation
-SIMID   = 'HD_study';  % Name of the simulation
+SIMID   = 'test_3D_marconi';  % Name of the simulation
+% SIMID   = 'HD_study';  % Name of the simulation
 % SIMID   = ['v3.0_P_',num2str(P),'_J_',num2str(J)];  % Name of the simulation
 NON_LIN = 1;   % activate non-linearity (is cancelled if KXEQ0 = 1)
 % INIT options
diff --git a/wk/plot_cosol_mat.m b/wk/plot_cosol_mat.m
index 876a233a..a6f583a9 100644
--- a/wk/plot_cosol_mat.m
+++ b/wk/plot_cosol_mat.m
@@ -1,29 +1,29 @@
-MAT = results.iCa;
-figure
-suptitle('DGGK,P=6,J=3');
-subplot(221)
-    imagesc(log(abs(MAT))); 
-    title('log abs')
-subplot(222)
-    imagesc(imag(MAT));  colormap(bluewhitered)
-    title('imag'); colorbar
-subplot(223)
-    imagesc(imag(MAT)>0);
-    title('imag$>$0');
-subplot(224)
-    imagesc(imag(MAT)<0);
-    title('imag$<$0');
-    
-    
+% MAT = results.iCa;
+% figure
+% suptitle('FCGK,P=6,J=3');
+% subplot(221)
+%     imagesc(log(abs(MAT))); 
+%     title('log abs')
+% subplot(222)
+%     imagesc(imag(MAT));  colormap(bluewhitered)
+%     title('imag'); colorbar
+% subplot(223)
+%     imagesc(imag(MAT)>0);
+%     title('imag$>$0');
+% subplot(224)
+%     imagesc(imag(MAT)<0);
+%     title('imag$<$0');
+%     
+%     
     %% SGGK
 P_ = 20; J_ = 10;
 mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5';
-SGGK_self = h5read(mat_file_name,'/00149/Caapj/Ciipj');
+SGGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj');
 SGGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
 SGGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
 
 figure
-MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=8');
+MAT = 1i*SGGK_self; suptitle('SGGK ii,P=20,J=10, k=0');
 % MAT = 1i*SGGK_ei; suptitle('SGGK ei,P=20,J=10');
 % MAT = 1i*SGGK_ie; suptitle('SGGK ie,P=20,J=10');
 subplot(221)
@@ -42,13 +42,38 @@ subplot(224)
         %% PAGK
 P_ = 20; J_ = 10;
 mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5';
-PAGK_self = h5read(mat_file_name,'/00149/Caapj/Ceepj');
+PAGK_self = h5read(mat_file_name,'/00000/Caapj/Ceepj');
 % PAGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
 % PAGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
 
 figure
-MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=8');
+MAT = 1i*PAGK_self; suptitle('PAGK ii,P=20,J=10, k=0');
+
+subplot(221)
+    imagesc(abs(MAT));
+    title('log abs')
+subplot(222)
+    imagesc(imag(MAT)); colormap(bluewhitered)
+    title('imag'); colorbar
+subplot(223)
+    imagesc(imag(MAT)>0);
+    title('imag$>$0');
+subplot(224)
+    imagesc(imag(MAT)<0);
+    title('imag$<$0');
 
+    
+    %% SGGK
+P_ = 6; J_ = 3;
+mat_file_name = '/home/ahoffman/HeLaZ/iCa/gk_coulomb_P_6_J_3_N_150_kpm_8.0.h5';
+FCGK_self = h5read(mat_file_name,'/00000/Caapj/Ciipj');
+FCGK_ei = h5read(mat_file_name,'/00000/Ceipj/CeipjF')+h5read(mat_file_name,'/00000/Ceipj/CeipjT');
+FCGK_ie = h5read(mat_file_name,'/00000/Ciepj/CiepjF')+h5read(mat_file_name,'/00000/Ciepj/CiepjT');
+
+figure
+MAT = 1i*FCGK_self; suptitle('FCGK ii,P=20,J=10, k=0');
+% MAT = 1i*SGGK_ei; suptitle('FCGK ei,P=20,J=10');
+% MAT = 1i*SGGK_ie; suptitle('FCGK ie,P=20,J=10');
 subplot(221)
     imagesc(abs(MAT));
     title('log abs')
@@ -61,3 +86,4 @@ subplot(223)
 subplot(224)
     imagesc(imag(MAT)<0);
     title('imag$<$0');
+    
\ No newline at end of file
-- 
GitLab