From 815be2d82380c23676ae61a1482275a44091ac51 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Mon, 5 Oct 2020 09:59:39 +0200
Subject: [PATCH] slight changes

---
 matlab/create_gif.m   |   7 ++-
 matlab/load_2D_data.m |   2 +-
 matlab/load_5D_data.m |   2 +-
 src/diagnose.F90      |   2 +-
 src/model_mod.F90     |   3 +-
 wk/analysis.m         | 141 ++++++++++++++++++++++++++++--------------
 wk/fort.90            |  40 ++++++------
 wk/setup.m            |  38 +++++++-----
 8 files changed, 146 insertions(+), 89 deletions(-)

diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 179a71e5..38bdbfed 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -23,9 +23,9 @@ fig  = figure;
     shading interp;
     hold on
   
-    n      = 0;
+    in      = 1;
     nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
-    for n = 1:numel(FIELD(1,1,:)) % time loop
+    for n = FRAMES % loop over selected frames
         pclr = pcolor(X,Y,FIELD(:,:,n)); shading interp; % frame plot
         set(pclr, 'edgecolor','none');
         title([FIELDNAME,', $t \approx$', sprintf('%.3d',ceil(T(n)))]);
@@ -35,7 +35,7 @@ fig  = figure;
         im = frame2im(frame); 
         [imind,cm] = rgb2ind(im,64); 
         % Write to the GIF File 
-        if n == 1 
+        if in == 1 
           imwrite(imind,cm,GIFNAME,'gif', 'Loopcount',inf); 
         else 
           imwrite(imind,cm,GIFNAME,'gif','WriteMode','append', 'DelayTime',DELAY); 
@@ -46,6 +46,7 @@ fig  = figure;
           nbytes = nbytes - 1;
         end
         nbytes = fprintf(2,'frame %d/%d',n,numel(FIELD(1,1,:)));
+        in = in + 1;
     end
     disp(' ')
     disp(['Gif saved @ : ',GIFNAME])
diff --git a/matlab/load_2D_data.m b/matlab/load_2D_data.m
index ddb42e1c..32c192fe 100644
--- a/matlab/load_2D_data.m
+++ b/matlab/load_2D_data.m
@@ -1,4 +1,4 @@
-function [ data, kr, kz, time ] = load_2D_data( filename, variablename )
+function [ data, kr, kz, time, dt ] = load_2D_data( filename, variablename )
 %LOAD_2D_DATA load a 2D variable stored in a hdf5 result file from HeLaZ
     time     = h5read(filename,'/data/var2d/time');
     kr       = h5read(filename,['/data/var2d/',variablename,'/coordkr']);
diff --git a/matlab/load_5D_data.m b/matlab/load_5D_data.m
index 34cb449f..767ca2b5 100644
--- a/matlab/load_5D_data.m
+++ b/matlab/load_5D_data.m
@@ -1,4 +1,4 @@
-function [ data, p, j, kr, kz, time ] = load_5D_data( filename, variablename )
+function [ data, p, j, kr, kz, time, dt ] = load_5D_data( filename, variablename )
 %LOAD_5D_DATA load a 5D variable stored in a hdf5 result file from HeLaZ
     time  = h5read(filename,'/data/var5d/time');
     p     = h5read(filename,['/data/var5d/',variablename,'/coordp']);
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 0c8a7dcb..7b217c7e 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -175,7 +175,7 @@ SUBROUTINE diagnose(kstep)
 
     ! Terminal info
     IF (MOD(cstep, INT(1.0/dt)) == 0) THEN
-     WRITE(*,"(F4.0,A,F4.0)") time,"/",tmax
+     WRITE(*,"(F5.0,A,F5.0)") time,"/",tmax
     ENDIF
 
      !                       2.1   0d history arrays
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 9d6f3ff9..58604cf5 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -6,6 +6,7 @@ MODULE model
   PRIVATE
 
   INTEGER,  PUBLIC, PROTECTED ::      CO =  0         ! Collision Operator
+  LOGICAL,  PUBLIC, PROTECTED ::      DK =  0         ! Drift kinetic model
   LOGICAL,  PUBLIC, PROTECTED :: NON_LIN =  .true.    ! To turn on non linear bracket term
   REAL(dp), PUBLIC, PROTECTED ::      mu =  0._dp     ! Hyperdiffusivity coefficient (for num. stability)
   REAL(dp), PUBLIC, PROTECTED ::      nu =  1._dp     ! Collision frequency
@@ -31,7 +32,7 @@ CONTAINS
     USE prec_const
     IMPLICIT NONE
 
-    NAMELIST /MODEL_PAR/ CO, NON_LIN, mu, nu, tau_e, tau_i, sigma_e, sigma_i, &
+    NAMELIST /MODEL_PAR/ CO, DK, NON_LIN, mu, nu, tau_e, tau_i, sigma_e, sigma_i, &
                          q_e, q_i, eta_n, eta_T, eta_B, lambdaD
 
     READ(lu_in,model_par)
diff --git a/wk/analysis.m b/wk/analysis.m
index 229dbd19..fbcac91a 100644
--- a/wk/analysis.m
+++ b/wk/analysis.m
@@ -2,14 +2,14 @@
 JOBNUM = 00;
 filename = [BASIC.SIMID,'_','%.2d.h5'];
 filename = sprintf(filename,JOBNUM); disp(['Analysing ',filename])
-[Nipj, p_, j_, kr, kz, Ts] = load_5D_data(filename, 'moments_i');
-Nepj                       = load_5D_data(filename, 'moments_e');
+[Nipj, p_, j_, kr, kz, Ts, dt] = load_5D_data(filename, 'moments_i');
+Nepj                           = load_5D_data(filename, 'moments_e');
 Ni      = squeeze(Nipj(1,1,:,:,:));
 Ne      = squeeze(Nepj(1,1,:,:,:));
 PH      = load_2D_data(filename, 'phi');
 Ts      = Ts';
 Ns      = numel(Ts);
-dt      = mean(diff(Ts));
+dt_samp = mean(diff(Ts));
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     Sipj    = load_5D_data(filename, 'Sipj');
     Sepj    = load_5D_data(filename, 'Sepj');
@@ -45,36 +45,54 @@ for it = 1:numel(PH(1,1,:))
 end
 
 %% Post processing
-phi_ST_r = zeros(Nr,Ns);   % Space-Time diagram of ES potential
-ne_ST_r  = zeros(Nr,Ns);   % Space-Time diagram of density
-ni_ST_r  = zeros(Nr,Ns);   % Space-Time diagram of density
-phi_ST_z = zeros(Nz,Ns);   % Space-Time diagram of ES potential
-ne_ST_z  = zeros(Nz,Ns);   % Space-Time diagram of density
-ni_ST_z  = zeros(Nz,Ns);   % Space-Time diagram of density
-phi_00 = zeros(1,Ns);    % Time evolution of ES potential at origin
-ne_00  = zeros(1,Ns);    % Time evolution of density at origin
-ni_00  = zeros(1,Ns);    % Time evolution of density at origin
-[~,ir0] = min(abs(r)); [~,iz0] = min(abs(z));
-Ne_11  = zeros(1,Ns);    % Time evolution of F density at 1,1
-Ni_11  = zeros(1,Ns);    % Time evolution of F density at 1,1
-[~,ikr1] = min(abs(kr-round(1/dkr)*dkr)); [~,ikz1] = min(abs(kz-round(1/dkz)*dkz));
-Sni_norm= zeros(1,Ns);   % Time evolution of the amp of density nonlin term
-Sne_norm= zeros(1,Ns);   % Time evolution of the amp of vorti. nonlin term
-E_pot  = zeros(1,Ns);    % Potential energy n^2
-E_kin  = zeros(1,Ns);    % Kinetic energy grad(phi)^2
-ExB    = zeros(1,Ns);    % ExB drift intensity \propto |\grad \phi|
-CFL    = zeros(1,Ns);    % CFL time step
+Ne_ST_kr = zeros(Nkr,Ns);   % Space-Time of max_kz Ne(k)
+Ne_ST_kz = zeros(Nkz,Ns);   % ''            max_kr Ne(k)
+ne_ST_r  = zeros(Nr,Ns);   % Space-Time of ne(z==0)
+ne_ST_z  = zeros(Nz,Ns);   % ''               r==0
+ne_00    = zeros(1,Ns);    % Time evolution of ne(r,z) at origin
+Ne_gm    = zeros(1,Ns);    % Time evolution of Ne(k) max gamma (max real)
+Ni_ST_kr = zeros(Nkr,Ns);   % same for ions
+Ni_ST_kz = zeros(Nkz,Ns);   % .
+ni_ST_r  = zeros(Nr,Ns);   % . 
+ni_ST_z  = zeros(Nz,Ns);   % .
+ni_00    = zeros(1,Ns);    % .
+Ni_gm    = zeros(1,Ns);    % .
+PH_ST_kr = zeros(Nkr,Ns);   % same for ES-potential
+PH_ST_kz = zeros(Nkz,Ns);   % . 
+phi_ST_r = zeros(Nr,Ns);   % .
+phi_ST_z = zeros(Nz,Ns);   % .
+phi_00   = zeros(1,Ns);    % .
+Sne_norm = zeros(1,Ns);    % Time evolution of the amp of e nonlin term
+Sni_norm = zeros(1,Ns);    % Time evolution of the amp of i nonlin term
+E_pot    = zeros(1,Ns);    % Potential energy n^2
+E_kin    = zeros(1,Ns);    % Kinetic energy grad(phi)^2
+ExB      = zeros(1,Ns);    % ExB drift intensity \propto |\grad \phi|
+CFL      = zeros(1,Ns);    % CFL time step
 Ddr = 1i*KR; Ddz = 1i*KZ; lapl   = Ddr.^2 + Ddz.^2; 
+[~,ir0]  = min(abs(r)); % index of r==0
+[~,iz0]  = min(abs(z)); % index of z==0
+[~,ikr1] = min(abs(kr-round(1/dkr)*dkr)); % index of kr==1 
+[~,ikz1] = min(abs(kz-round(1/dkz)*dkz)); % index of kz==1
 for it = 1:numel(PH(1,1,:))
     NE_ = Ne(:,:,it); NN_ = Ni(:,:,it); PH_ = PH(:,:,it);
-    phi_ST_r(:,it) = phi(:,iz0,it); phi_ST_z(:,it) = phi(ir0,:,it);
+
     ne_ST_r  (:,it)= ne(:,iz0,it);  ne_ST_z  (:,it)= ne(ir0,:,it);
+    Ne_ST_kr (:,it)= max(real(Ne(:,:,it)),[],2); 
+    Ne_ST_kz (:,it)= max(real(Ne(:,:,it)),[],1);
+    ne_00(it)      = ne(ir0,iz0,it);
+    Ne_gm(it)      = max(max(real(Ne(:,:,it))));
+    
     ni_ST_r  (:,it)= ni(:,iz0,it);  ni_ST_z  (:,it)= ni(ir0,:,it);
-    phi_00(it)   = phi(ir0,iz0,it);
-    ne_00(it)    = ne(ir0,iz0,it);
-    ni_00(it)    = ni(ir0,iz0,it);
-    Ne_11(it)    = Ne(ikr1,ikz1,it);
-    Ni_11(it)    = Ni(ikr1,ikz1,it);
+    Ni_ST_kr (:,it)= max(real(Ni(:,:,it)),[],2); 
+    Ni_ST_kz (:,it)= max(real(Ni(:,:,it)),[],1);
+    ni_00(it)      = ni(ir0,iz0,it);
+    Ni_gm(it)      = max(max(real(Ni(:,:,it))));
+    
+    phi_ST_r(:,it) = phi(:,iz0,it); phi_ST_z(:,it) = phi(ir0,:,it);
+    PH_ST_kr(:,it) = max(real(PH(:,:,it)),[],2); 
+    PH_ST_kz(:,it) = max(real(PH(:,:,it)),[],1);
+    phi_00(it)     = phi(ir0,iz0,it);
+    
 if strcmp(OUTPUTS.write_non_lin,'.true.')
     Sni_norm(it) = sum(sum(abs(SNi(:,:,it))));
     Sne_norm(it) = sum(sum(abs(SNe(:,:,it))));
@@ -94,7 +112,7 @@ fig = figure; FIGNAME = ['t_evolutions',sprintf('_%.2d',JOBNUM)];
         semilogy(Ts,abs(ni_00),'-','DisplayName','$n_i^{00}$');
         grid on; xlabel('$t$'); ylabel('$|n_a(x=0,y=0)|$');
     subplot(222)
-        semilogy(Ts,abs(Ni_11),'-','DisplayName','$\phi$')
+        semilogy(Ts,abs(Ni_gm),'-','DisplayName','$\phi$')
         grid on; xlabel('$t$'); ylabel('$|\tilde n(k_r\approx 1,k_z\approx 1)|$');
     subplot(223)
         semilogy(Ts,E_kin+E_pot,'-','DisplayName','$\sum|ik\tilde\phi_i|^2+\sum|\tilde n_i|^2$')
@@ -110,12 +128,12 @@ end
 FMT = '.fig'; save_figure
 
 %% Spectra energy
-fig = figure; FIGNAME = ['Energy_kin_KZ',sprintf('_%.2d',JOBNUM)];
-    semilogy(kr(floor(end/2)+1:end),E_kin_KR(floor(end/2)+1:end),'o','DisplayName','$\sum_y\langle|ik\tilde\phi_i|^2\rangle_t$')
-    hold on;
-    loglog(kz(floor(end/2)+1:end),E_kin_KZ(floor(end/2)+1:end),'o','DisplayName','$\sum_x\langle|ik\tilde\phi_i|^2\rangle_t$')
-    grid on; xlabel('$k$');  legend('show');
-FMT = '.fig'; save_figure
+% fig = figure; FIGNAME = ['Energy_kin_KZ',sprintf('_%.2d',JOBNUM)];
+%     semilogy(kr(floor(end/2)+1:end),E_kin_KR(floor(end/2)+1:end),'o','DisplayName','$\sum_y\langle|ik\tilde\phi_i|^2\rangle_t$')
+%     hold on;
+%     loglog(kz(floor(end/2)+1:end),E_kin_KZ(floor(end/2)+1:end),'o','DisplayName','$\sum_x\langle|ik\tilde\phi_i|^2\rangle_t$')
+%     grid on; xlabel('$k$');  legend('show');
+% FMT = '.fig'; save_figure
 
 %% CFL condition
 fig = figure; FIGNAME = ['CFL',sprintf('_%.2d',JOBNUM)];
@@ -155,17 +173,41 @@ subplot(223)% density
     xlabel('$z\,(r=0)$'); ylabel('$t$'); title('$\phi$');
 FMT = '.fig'; save_figure
 
-%% phi
-fig = figure; FIGNAME = ['phi_ST',sprintf('_%.2d',JOBNUM)];
-    [TY,TX] = meshgrid(Ts,z);
-    pclr = pcolor(TX,TY,(plt(phi_ST_r))); set(pclr, 'edgecolor','none'); colorbar;
-    xlabel('$x\,(y=0)$'); ylabel('$t$'); title('$\phi$');
+%% Space-Time diagrams for max_kz(Real)
+plt = @(x) (x);
+fig = figure; FIGNAME = ['kr_space_time_diag',sprintf('_%.2d',JOBNUM)];
+    [TY,TX] = meshgrid(Ts,kr);
+subplot(221)% density
+    pclr = pcolor(TX,TY,(plt(Ne_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(N_e^{00}))$');
+subplot(222)% density
+    pclr = pcolor(TX,TY,(plt(Ni_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(N_i^{00}))$');
+subplot(223)% density
+    pclr = pcolor(TX,TY,(plt(PH_ST_kr))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kr$'); ylabel('$t$'); title('$\max_{k_z}(\textrm{Re}(\tilde\phi$))');
 FMT = '.fig'; save_figure
 
+%% Space-Time diagrams at max_kr(Real)
+plt = @(x) (x);
+fig = figure; FIGNAME = ['kz_space_time_diag',sprintf('_%.2d',JOBNUM)];
+    [TY,TX] = meshgrid(Ts,kz);
+subplot(221)% density
+    pclr = pcolor(TX,TY,(plt(Ne_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(N_e^{00}))$');
+subplot(222)% density
+    pclr = pcolor(TX,TY,(plt(Ni_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(N_i^{00}))$');
+subplot(223)% density
+    pclr = pcolor(TX,TY,(plt(PH_ST_kz))); set(pclr, 'edgecolor','none'); colorbar;
+    xlabel('$kz$'); ylabel('$t$'); title('$\max_{k_r}(\textrm{Re}(\tilde\phi))$');
+FMT = '.fig'; save_figure
+
+
 if 0
 %% Show frame
-it = min(1,numel(Ts));
-fig = figure; FIGNAME = ['frame',sprintf('_%.2d',SID)];
+it = min(50,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;
         xlabel('$k_r$'); ylabel('$k_z$'); title(sprintf('t=%.3d',Ts(it))); legend('$\hat\phi$');
@@ -183,19 +225,24 @@ end
 FMT = '.fig'; save_figure
 end
 %%
-DELAY = 0.1; skip_ = 2;
-if 0
+DELAY = 0.07; skip_ = 1;
+FRAMES = 200:skip_:numel(Ts);
+if 1
 %% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Density electron
 GIFNAME = ['ne',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_e^{00}$';
-FIELD = real(ne(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+FIELD = real(ne); X = XX; Y = YY; T = Ts;
 create_gif
 %% Density ion
 GIFNAME = ['ni',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$n_i^{00}$';
-FIELD = real(ni(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+FIELD = real(ni); X = XX; Y = YY; T = Ts;
 create_gif
 %% Phi
 GIFNAME = ['phi',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$\phi$';
-FIELD = real(phi(:,:,1:skip_:end)); X = XX; Y = YY; T = Ts(1:skip_:end);
+FIELD = real(phi); X = XX; Y = YY; T = Ts;
+create_gif
+%% Density electron frequency
+GIFNAME = ['Ni',sprintf('_%.2d',JOBNUM)]; FIELDNAME = '$N_i^{00}$';
+FIELD = real(Ni); X = KR; Y = KZ; T = Ts;
 create_gif
 end
\ No newline at end of file
diff --git a/wk/fort.90 b/wk/fort.90
index c11d2956..7001202e 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,18 +1,19 @@
 &BASIC
   nrun   = 100000000
-  dt     = 0.001
-  tmax   = 200
+  dt     = 0.0001
+  tmax   = 100
   RESTART = .true.
 /
 &GRID
-  pmaxe =1
-  jmaxe = 0
-  pmaxi = 1
-  jmaxi = 0
+  pmaxe =5
+  jmaxe = 2
+  pmaxi = 5
+  jmaxi = 2
   Nr   = 64
-  Lr = 20
+  Lr = 10
   Nz   = 64
-  Lz = 20
+  Lz = 10
+  kpar = 0
 /
 &OUTPUT_PAR
   nsave_0d = 0
@@ -22,18 +23,19 @@
   write_Ni00    = .true.
   write_moments = .true.
   write_phi     = .true.
-  write_non_lin     = .true.
+  write_non_lin     = .false.
   write_doubleprecision = .true.
-  resfile0      = 'SwoK_32x64_L_20_P_1_J_0_nB_0.5_nN_1_mu_1e-04_'
-  rstfile0      = '../checkpoint/cp_SwoK_32x64_L_20_P_1_J_0_nB_0.5_nN_1_mu_1e-04_'
+  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_'
   job2load      = 0
 /
 &MODEL_PAR
   ! Collisionality
-  CO      = 0
-  NON_LIN = .true.
-  mu      = 0.0001
-  nu      = 0
+  CO      = -2
+  DK      = .false.
+  NON_LIN = .false.
+  mu      = 0.01
+  nu      = 0.01
   tau_e   = 1
   tau_i   = 1
   sigma_e = 0.023338
@@ -42,7 +44,7 @@
   q_i     = 1
   eta_n   = 1
   eta_T   = 0
-  eta_B   = 0.5
+  eta_B   = 0.1
   lambdaD = 0
 /
 &INITIAL_CON
@@ -50,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_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'
+  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'
 /
 &TIME_INTEGRATION_PAR
   numerical_scheme='RK4'
diff --git a/wk/setup.m b/wk/setup.m
index c1c45262..022ec000 100644
--- a/wk/setup.m
+++ b/wk/setup.m
@@ -5,29 +5,32 @@ default_plots_options
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.0;    % Collision frequency
+NU      = 1e-2;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.5;    % Magnetic gradient
+ETAB    = 0.1;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
-MU      = 1e-4;   % Hyper diffusivity coefficient
+MU      = 1e-2;   % Hyper diffusivity coefficient
+LAMBDAD = 0.0;
 %% GRID PARAMETERS
-N       = 64;     % Frequency gridpoints
-L       = 20;     % Size of the squared frequency domain
-PMAXE   = 01;     % Highest electron Hermite polynomial degree
-JMAXE   = 00;     % Highest ''       Laguerre ''
-PMAXI   = 01;     % Highest ion      Hermite polynomial degree
-JMAXI   = 00;     % Highest ''       Laguerre ''
+N       = 64;     % Frequency gridpoints (Nr = 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 ''
+KPAR    = 0.0;    % Parellel wave vector component
 %% TIME PARAMETERS 
-TMAX    = 200.0;    % Maximal time unit
-DT      = 1e-3;   % Time step
-SPS     = 1;     % Sampling per time unit
+TMAX    = 100.0;    % Maximal time unit
+DT      = 1e-4;   % Time step
+SPS     = 10;     % Sampling per time unit
 RESTART = 1;      % To restart from last checkpoint
 JOB2LOAD= 0;
 %% OPTIONS
-SIMID   = 'SwoK';  % Name of the simulation
-NON_LIN = 1; % activate non-linearity 
-CO      = 0;
+SIMID   = 'gvskr';  % 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)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -37,6 +40,7 @@ params   = ['P_',num2str(PMAXE),'_J_',num2str(JMAXE),...
     '_nB_',num2str(ETAB),'_nN_',num2str(ETAN),'_mu_%0.0e_'];
 params   = sprintf(params,MU);
 if ~NON_LIN; params = ['lin_',params]; end;
+if  DK;      params = ['DK_',params]; end;
 resolution = ['_',num2str(N/2),'x',num2str(N),'_'];
 gridname   = ['L_',num2str(L),'_'];
 BASIC.SIMID = [SIMID,resolution,gridname,params];
@@ -66,8 +70,10 @@ GRID.Nr    = N; % r grid resolution
 GRID.Lr    = L; % r length
 GRID.Nz    = N; % z ''
 GRID.Lz    = L; % z ''
+GRID.kpar  = KPAR;
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+if DK;      MODEL.DK      = '.true.'; else; MODEL.DK      = '.false.';end;
 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
@@ -84,7 +90,7 @@ MODEL.q_i     = 1.0;
 MODEL.eta_n   = ETAN;        % source term kappa for HW
 MODEL.eta_T   = ETAT;        % Temperature
 MODEL.eta_B   = ETAB;        % Magnetic
-MODEL.lambdaD = 0.0;
+MODEL.lambdaD = LAMBDAD;
 % Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
 INITIAL.only_Na00         = '.false.';
-- 
GitLab