From 06b78cf31a8778bfc19bbc50518b1c7005ac4b42 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 5 Oct 2020 10:20:03 +0200
Subject: [PATCH] put kz > 0 instead of kr > 0

---
 matlab/create_gif.m |  3 ++-
 src/grid_mod.F90    | 17 +++++++++++------
 wk/analysis.m       | 22 +++++++++++++---------
 wk/fort.90          | 28 ++++++++++++++--------------
 wk/setup.m          | 16 ++++++++--------
 5 files changed, 48 insertions(+), 38 deletions(-)

diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 38bdbfed..5b00109e 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -21,10 +21,11 @@ fig  = figure;
     colorbar
     axis tight manual % this ensures that getframe() returns a consistent size
     shading interp;
+    xlabel(XNAME); ylabel(YNAME);
     hold on
   
     in      = 1;
-    nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+    nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,:)));
     for n = FRAMES % loop over selected frames
         pclr = pcolor(X,Y,FIELD(:,:,n)); shading interp; % frame plot
         set(pclr, 'edgecolor','none');
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 0987ba28..6fee2179 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -128,7 +128,8 @@ CONTAINS
     IMPLICIT NONE
     INTEGER :: ikr
 
-    Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
+    ! Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
+    Nkr = Nr
     ! Start and END indices of grid
     ikrs = 1
     ikre = Nkr
@@ -138,11 +139,13 @@ CONTAINS
     ! Discretized kr positions ordered as dk*(0 1 2)
     ALLOCATE(krarray(ikrs:ikre))
     DO ikr = ikrs,ikre
-      krarray(ikr) = REAL(ikr-1,dp) * deltakr
+      ! krarray(ikr) = REAL(ikr-1,dp) * deltakr
+      krarray(ikr) = deltakr*(MODULO(ikr-1,Nkr/2)-Nkr/2*FLOOR(2.*real(ikr-1)/real(Nkr)))
       if (krarray(ikr) .EQ. 0) THEN
         ikr_0 = ikr
       ENDIF
     END DO
+    krarray(Nr/2+1) = -krarray(Nr/2+1)
 
     ! Orszag 2/3 filter
     two_third_krmax = 2._dp/3._dp*deltakr*Nkr
@@ -160,22 +163,24 @@ CONTAINS
     USE prec_const
     IMPLICIT NONE
 
-    Nkz = Nz;
+    ! Nkz = Nz;
+    Nkz = Nz/2 + 1; ! Defined only on positive kz since fields are real
     ! Start and END indices of grid
     ikzs = 1
     ikze = Nkz
     ! Grid spacings
-    deltakz = 2._dp*PI/Nkz/deltaz
+    deltakz = 2._dp*PI/Nz/deltaz
 
     ! Discretized kz positions ordered as dk*(0 1 2 -3 -2 -1)
     ALLOCATE(kzarray(ikzs:ikze))
     DO ikz = ikzs,ikze
-      kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
+      ! kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
+      kzarray(ikz) = REAL(ikz-1,dp) * deltakz
       if (kzarray(ikz) .EQ. 0) THEN
         ikz_0 = ikz
       ENDIF
     END DO
-    kzarray(Nz/2+1) = -kzarray(Nz/2+1)
+    ! kzarray(Nz/2+1) = -kzarray(Nz/2+1)
 
     ! Orszag 2/3 filter
     two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
diff --git a/wk/analysis.m b/wk/analysis.m
index fbcac91a..84232c45 100644
--- a/wk/analysis.m
+++ b/wk/analysis.m
@@ -206,7 +206,7 @@ FMT = '.fig'; save_figure
 
 if 0
 %% Show frame
-it = min(50,numel(Ts));
+it = min(70,numel(Ts));
 fig = figure; FIGNAME = ['frame',sprintf('_%.2d',JOBNUM)];
     subplot(221); plt = @(x) fftshift((real(x)));
         pclr = pcolor(fftshift(KR),fftshift(KZ),plt(PH(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
@@ -225,24 +225,28 @@ end
 FMT = '.fig'; save_figure
 end
 %%
-DELAY = 0.07; skip_ = 1;
-FRAMES = 200:skip_:numel(Ts);
-if 1
+DELAY = 0.07; skip_ = 10;
+FRAMES = 1:skip_:numel(Ts);
+if 0
 %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Density electron
-GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_e^{00}$';
+GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)];
 FIELD = real(ne); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$n_e^{00}$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Density ion
-GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_i^{00}$';
+GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; 
 FIELD = real(ni); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$n_i^{00}$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Phi
-GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$\phi$';
+GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)]; 
 FIELD = real(phi); X = XX; Y = YY; T = Ts;
+FIELDNAME = '$\phi$'; XNAME = '$r$'; YNAME = '$z$';
 create_gif
 %% Density electron frequency
-GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$N_i^{00}$';
-FIELD = real(Ni); X = KR; Y = KZ; T = Ts;
+GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)]; 
+FIELD = fftshift(real(Ni)); X = fftshift(KR); Y = fftshift(KZ); T = Ts;
+FIELDNAME = '$N_i^{00}$'; XNAME = '$k_r$'; YNAME = '$k_z$';
 create_gif
 end
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index 7001202e..87909eb5 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,14 +1,14 @@
 &BASIC
   nrun   = 100000000
-  dt     = 0.0001
-  tmax   = 100
-  RESTART = .true.
+  dt     = 0.01
+  tmax   = 10
+  RESTART = .false.
 /
 &GRID
-  pmaxe =5
-  jmaxe = 2
-  pmaxi = 5
-  jmaxi = 2
+  pmaxe =1
+  jmaxe = 0
+  pmaxi = 1
+  jmaxi = 0
   Nr   = 64
   Lr = 10
   Nz   = 64
@@ -18,15 +18,15 @@
 &OUTPUT_PAR
   nsave_0d = 0
   nsave_1d = 0
-  nsave_2d = 1000
-  nsave_5d = 1000
+  nsave_2d = 10
+  nsave_5d = 10
   write_Ni00    = .true.
   write_moments = .true.
   write_phi     = .true.
   write_non_lin     = .false.
   write_doubleprecision = .true.
-  resfile0      = 'gvskr_32x64_L_10_lin_P_5_J_2_nB_0.1_nN_1_mu_1e-02_'
-  rstfile0      = '../checkpoint/cp_gvskr_32x64_L_10_lin_P_5_J_2_nB_0.1_nN_1_mu_1e-02_'
+  resfile0      = 'test_kzpos_32x64_L_10_lin_P_1_J_0_nB_0.1_nN_1_mu_1e-02_'
+  rstfile0      = '../checkpoint/cp_test_kzpos_32x64_L_10_lin_P_1_J_0_nB_0.1_nN_1_mu_1e-02_'
   job2load      = 0
 /
 &MODEL_PAR
@@ -52,9 +52,9 @@
   initback_moments  =0.0001
   initnoise_moments =5e-05
   iseed             =42
-  selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_5_Jmaxe_2_Pmaxi_5_Jmaxi_2_pamaxx_10.h5'
-  eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_5_Jmaxe_2_Pmaxi_5_Jmaxi_2_pamaxx_10_tau_1.0000_mu_0.0233.h5'
-  iemat_file        ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_5_Jmaxe_2_Pmaxi_5_Jmaxi_2_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+  selfmat_file      ='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10.h5'
+  eimat_file        ='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10_tau_1.0000_mu_0.0233.h5'
+  iemat_file        ='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_1_Jmaxe_0_Pmaxi_1_Jmaxi_0_pamaxx_10_tau_1.0000_mu_0.0233.h5'
 /
 &TIME_INTEGRATION_PAR
   numerical_scheme='RK4'
diff --git a/wk/setup.m b/wk/setup.m
index a89dd6e6..80bdb5c7 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -15,19 +15,19 @@ LAMBDAD = 0.0;
 %% GRID PARAMETERS
 N       = 64;     % Frequency gridpoints (Nkr = N/2)
 L       = 10;      % Size of the squared frequency domain
-PMAXE   = 05;     % Highest electron Hermite polynomial degree
-JMAXE   = 02;     % Highest ''       Laguerre ''
-PMAXI   = 05;     % Highest ion      Hermite polynomial degree
-JMAXI   = 02;     % Highest ''       Laguerre ''
+PMAXE   = 01;     % Highest electron Hermite polynomial degree
+JMAXE   = 00;     % Highest ''       Laguerre ''
+PMAXI   = 01;     % Highest ion      Hermite polynomial degree
+JMAXI   = 00;     % Highest ''       Laguerre ''
 KPAR    = 0.0;    % Parellel wave vector component
 %% TIME PARAMETERS 
-TMAX    = 100.0;    % Maximal time unit
-DT      = 1e-4;   % Time step
+TMAX    = 10.0;    % Maximal time unit
+DT      = 1e-2;   % Time step
 SPS     = 10;     % Sampling per time unit
-RESTART = 1;      % To restart from last checkpoint
+RESTART = 0;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = 'gvskr';  % Name of the simulation
+SIMID   = 'test_kzpos';  % Name of the simulation
 NON_LIN = 0;   % activate non-linearity (is cancelled if KREQ0 = 1)
 CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
 DK      = 0;   % Drift kinetic model (put every kernel to 1)
-- 
GitLab