From 27afc84e76533c8fa63df8c8d628ad1fb0e2b53f Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 3 Aug 2021 14:26:05 +0200
Subject: [PATCH] added shear and inverse aspect ratio as input parameters

---
 matlab/setup.m        |  6 +++--
 matlab/write_fort90.m | 52 ++++++++++++++++++++++---------------------
 src/grid_mod.F90      | 36 ++++++++++++++++--------------
 wk/analysis_3D.m      |  6 ++---
 wk/local_run.m        | 25 ++++++++++++---------
 wk/marconi_run.m      | 12 +++++-----
 6 files changed, 73 insertions(+), 64 deletions(-)

diff --git a/matlab/setup.m b/matlab/setup.m
index dd614b74..0f6ded1f 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -9,8 +9,10 @@ GRID.Nx    = N; % x grid resolution
 GRID.Lx    = L; % x length
 GRID.Ny    = N * (1-KXEQ0) + KXEQ0; % y ''
 GRID.Ly    = L * (1-KXEQ0); % y ''
-GRID.Nz    = Nz; % z resolution
-GRID.q0    = q0; % q factor
+GRID.Nz    = Nz;    % z resolution
+GRID.q0    = q0;    % q factor
+GRID.shear = shear; % shear
+GRID.eps   = eps;   % inverse aspect ratio
 
 % Model parameters
 MODEL.CO      = CO;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 366ee163..311229af 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -7,21 +7,23 @@ fprintf(fid,'&BASIC\n');
 fprintf(fid,['  nrun   = ', num2str(BASIC.nrun),'\n']);
 fprintf(fid,['  dt     = ', num2str(BASIC.dt),'\n']);
 fprintf(fid,['  tmax   = ', num2str(BASIC.tmax),'\n']);
-fprintf(fid,['  RESTART   = ', BASIC.RESTART,'\n']);
+fprintf(fid,['  RESTART    = ', BASIC.RESTART,'\n']);
 fprintf(fid,['  maxruntime = ', num2str(BASIC.maxruntime),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&GRID\n');
-fprintf(fid,['  pmaxe =', num2str(GRID.pmaxe),'\n']);
-fprintf(fid,['  jmaxe = ', num2str(GRID.jmaxe),'\n']);
-fprintf(fid,['  pmaxi = ', num2str(GRID.pmaxi),'\n']);
-fprintf(fid,['  jmaxi = ', num2str(GRID.jmaxi),'\n']);
-fprintf(fid,['  Nx   = ', num2str(GRID.Nx),'\n']);
-fprintf(fid,['  Lx = ', num2str(GRID.Lx),'\n']);
-fprintf(fid,['  Ny   = ', num2str(GRID.Ny),'\n']);
-fprintf(fid,['  Ly = ', num2str(GRID.Ly),'\n']);
-fprintf(fid,['  Nz   = ', num2str(GRID.Nz),'\n']);
-fprintf(fid,['  q0  = ', num2str(GRID.q0),'\n']);
+fprintf(fid,['  pmaxe  =', num2str(GRID.pmaxe),'\n']);
+fprintf(fid,['  jmaxe  = ', num2str(GRID.jmaxe),'\n']);
+fprintf(fid,['  pmaxi  = ', num2str(GRID.pmaxi),'\n']);
+fprintf(fid,['  jmaxi  = ', num2str(GRID.jmaxi),'\n']);
+fprintf(fid,['  Nx     = ', num2str(GRID.Nx),'\n']);
+fprintf(fid,['  Lx     = ', num2str(GRID.Lx),'\n']);
+fprintf(fid,['  Ny     = ', num2str(GRID.Ny),'\n']);
+fprintf(fid,['  Ly     = ', num2str(GRID.Ly),'\n']);
+fprintf(fid,['  Nz     = ', num2str(GRID.Nz),'\n']);
+fprintf(fid,['  q0     = ', num2str(GRID.q0),'\n']);
+fprintf(fid,['  shear  = ', num2str(GRID.q0),'\n']);
+fprintf(fid,['  eps    = ', num2str(GRID.q0),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&OUTPUT_PAR\n');
@@ -34,11 +36,11 @@ fprintf(fid,['  nsave_cp = ', num2str(OUTPUTS.nsave_cp),'\n']);
 fprintf(fid,['  write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
 fprintf(fid,['  write_gamma = ', OUTPUTS.write_gamma,'\n']);
 fprintf(fid,['  write_phi   = ', OUTPUTS.write_phi,'\n']);
-fprintf(fid,['  write_Na00 = ', OUTPUTS.write_Na00,'\n']);
-fprintf(fid,['  write_Napj = ', OUTPUTS.write_Napj,'\n']);
-fprintf(fid,['  write_Sapj = ', OUTPUTS.write_Sapj,'\n']);
-fprintf(fid,['  write_dens = ', OUTPUTS.write_dens,'\n']);
-fprintf(fid,['  write_temp = ', OUTPUTS.write_temp,'\n']);
+fprintf(fid,['  write_Na00  = ', OUTPUTS.write_Na00,'\n']);
+fprintf(fid,['  write_Napj  = ', OUTPUTS.write_Napj,'\n']);
+fprintf(fid,['  write_Sapj  = ', OUTPUTS.write_Sapj,'\n']);
+fprintf(fid,['  write_dens  = ', OUTPUTS.write_dens,'\n']);
+fprintf(fid,['  write_temp  = ', OUTPUTS.write_temp,'\n']);
 fprintf(fid,['  resfile0      = ', OUTPUTS.resfile0,'\n']);
 fprintf(fid,['  rstfile0      = ', OUTPUTS.rstfile0,'\n']);
 fprintf(fid,['  job2load      = ', num2str(OUTPUTS.job2load),'\n']);
@@ -67,18 +69,18 @@ fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&INITIAL_CON\n');
-fprintf(fid,['  INIT_NOISY_PHI    =', INITIAL.init_noisy_phi,'\n']);
-fprintf(fid,['  INIT_ZF       =', num2str(INITIAL.INIT_ZF),'\n']);
-fprintf(fid,['  WIPE_TURB     =', INITIAL.wipe_turb,'\n']);
-fprintf(fid,['  INIT_BLOB     =', INITIAL.init_blob,'\n']);
-fprintf(fid,['  init_background  =', num2str(INITIAL.init_background),'\n']);
-fprintf(fid,['  init_noiselvl =', num2str(INITIAL.init_noiselvl),'\n']);
-fprintf(fid,['  iseed             =', num2str(INITIAL.iseed),'\n']);
-fprintf(fid,['  mat_file        =', INITIAL.mat_file,'\n']);
+fprintf(fid,['  INIT_NOISY_PHI    = ', INITIAL.init_noisy_phi,'\n']);
+fprintf(fid,['  INIT_ZF       = ', num2str(INITIAL.INIT_ZF),'\n']);
+fprintf(fid,['  WIPE_TURB     = ', INITIAL.wipe_turb,'\n']);
+fprintf(fid,['  INIT_BLOB     = ', INITIAL.init_blob,'\n']);
+fprintf(fid,['  init_background  = ', num2str(INITIAL.init_background),'\n']);
+fprintf(fid,['  init_noiselvl = ', num2str(INITIAL.init_noiselvl),'\n']);
+fprintf(fid,['  iseed         = ', num2str(INITIAL.iseed),'\n']);
+fprintf(fid,['  mat_file      = ', INITIAL.mat_file,'\n']);
 fprintf(fid,'/\n');
 
 fprintf(fid,'&TIME_INTEGRATION_PAR\n');
-fprintf(fid,['  numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
+fprintf(fid,['  numerical_scheme = ', TIME_INTEGRATION.numerical_scheme,'\n']);
 fprintf(fid,'/');
 
 fclose(fid);
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index f925d1bc..42d612c1 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -82,6 +82,19 @@ MODULE grid
 
 CONTAINS
 
+
+  SUBROUTINE grid_readinputs
+    ! Read the input parameters
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: lu_in   = 90              ! File duplicated from STDIN
+
+    NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
+                    Nx,  Lx,  Ny,  Ly, Nz, q0, shear, eps
+    READ(lu_in,grid)
+
+  END SUBROUTINE grid_readinputs
+
   SUBROUTINE init_1Dgrid_distr
 
     ! write(*,*) Nx
@@ -288,19 +301,6 @@ CONTAINS
     if(my_id.EQ.0) write(*,*) '#parallel planes = ', Nz
   END SUBROUTINE set_zgrid
 
-  SUBROUTINE grid_readinputs
-    ! Read the input parameters
-    USE prec_const
-    IMPLICIT NONE
-    INTEGER :: lu_in   = 90              ! File duplicated from STDIN
-
-    NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
-                    Nx,  Lx,  Ny,  Ly, Nz, q0
-    READ(lu_in,grid)
-
-  END SUBROUTINE grid_readinputs
-
-
   SUBROUTINE grid_outputinputs(fidres, str)
     ! Write the input parameters to the results_xx.h5 file
 
@@ -321,10 +321,12 @@ CONTAINS
     CALL attach(fidres, TRIM(str),   "Ly",   Ly)
     CALL attach(fidres, TRIM(str),   "Nz",   Nz)
     CALL attach(fidres, TRIM(str),   "q0",   q0)
-    CALL attach(fidres, TRIM(str),   "Nkx",   Nkx)
-    CALL attach(fidres, TRIM(str),   "Lkx",   Lkx)
-    CALL attach(fidres, TRIM(str),   "Nky",   Nky)
-    CALL attach(fidres, TRIM(str),   "Lky",   Lky)
+    CALL attach(fidres, TRIM(str),"shear",shear)
+    CALL attach(fidres, TRIM(str),  "eps",  eps)
+    CALL attach(fidres, TRIM(str),  "Nkx",  Nkx)
+    CALL attach(fidres, TRIM(str),  "Lkx",  Lkx)
+    CALL attach(fidres, TRIM(str),  "Nky",  Nky)
+    CALL attach(fidres, TRIM(str),  "Lky",  Lky)
   END SUBROUTINE grid_outputinputs
 
   FUNCTION bare(p_,j_)
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index e146bc66..5300342d 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -7,10 +7,8 @@ outfile ='';
 outfile ='';
 outfile ='';
 outfile ='';
-outfile ='';
-outfile ='';
-outfile ='HD_study/150x75_L_100_P_4_J_2_eta_0.6_nu_1e-02_DGGK_mu_3e-03';
-% outfile ='HD_study/300x150_L_100_P_2_J_1_eta_0.6_nu_5e-01_DGGK_mu_2e-03';
+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';
     BASIC.RESDIR      = ['../results/',outfile,'/'];
     BASIC.MISCDIR     = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
     CMD = ['cp ', BASIC.RESDIR,'outputs* ',BASIC.MISCDIR]; disp(CMD);
diff --git a/wk/local_run.m b/wk/local_run.m
index 1574521b..c588c465 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -4,18 +4,18 @@ addpath(genpath('../matlab')) % ... add
 CLUSTER.TIME  = '99:00:00'; % allocation time hh:mm:ss
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
-NU      = 0.5;   % Collision frequency
+NU      = 0.1;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
-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)
-q0      = 1.0;    % q factor ()
+NU_HYP  = 0.1;
+%% GRID AND GEOMETRY PARAMETERS
+N       = 100;     % Frequency gridpoints (Nkx = N/2)
+L       = 60;     % Size of the squared frequency domain
+Nz      = 10;      % number of perpendicular planes (parallel grid)
+q0      = 1.0;    % safety factor
+shear   = 0.0;    % magnetic shear
+eps     = 0.0;    % inverse aspect ratio
 P       = 2;
 J       = 1;
-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 PARAMETERS
 TMAX    = 200;  % Maximal time unit
 DT      = 1e-2;   % Time step
@@ -32,12 +32,13 @@ JOB2LOAD= 0;
 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   = 'HD_study';  % Name of the simulation
+% SIMID   = 'HD_study';  % Name of the simulation
+SIMID   = 'test_3D';  % 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
 INIT_ZF = 0; ZF_AMP = 0.0;
-INIT_BLOB = 1; WIPE_TURB = 1;
+INIT_BLOB = 0; WIPE_TURB = 0;
 %% OUTPUTS
 W_DOUBLE = 0;
 W_GAMMA  = 1;
@@ -68,6 +69,8 @@ TAU     = 1.0;    % e/i temperature ratio
 ETAT    = 0.0;    % Temperature gradient
 ETAB    = 1.0;    % Magnetic gradient (1.0 to set R=LB)
 INIT_PHI= 1;   % Start simulation with a noisy phi and moments
+MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 %% Setup and file management
 setup
 system('rm fort.90');
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 3b8f42ed..236dc473 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -2,7 +2,7 @@ clear all;
 addpath(genpath('../matlab')) % ... add
 SUBMIT = 1; % To submit the job automatically
 % EXECNAME = 'helaz_dbg';
-  EXECNAME = 'helaz_3.0';
+  EXECNAME = 'helaz_3.1';
 for ETAN = [1/0.6]
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
@@ -20,16 +20,16 @@ NP_KX         = 24;         % MPI processes along kx
 %% PHYSICAL PARAMETERS
 NU      = 0.05;   % Collision frequency
 ETAN    = 1.0/0.6;    % Density gradient drive (R/Ln)
-NU_HYP  = 1.0;
+NU_HYP  = 10.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)
 q0      = 1.0;    % q factor ()
+shear   = 0.0;    % magnetic shear
+eps     = 0.0;    % inverse aspect ratio
 P       = 2;
 J       = 1;
-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 PARAMETERS
 TMAX    = 200;  % Maximal time unit
 DT      = 1e-2;   % Time step
@@ -43,7 +43,7 @@ JOB2LOAD= 0;
 %% OPTIONS AND NAMING
 % Collision operator
 % (0 : L.Bernstein, 1 : Dougherty, 2: Sugama, 3 : Pitch angle ; +/- for GK/DK)
-CO      = 1;
+CO      = 2;
 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
@@ -83,6 +83,8 @@ TAU     = 1.0;    % e/i temperature ratio
 ETAT    = 0.0;    % Temperature gradient
 ETAB    = 1.0;    % Magnetic gradient (1.0 to set R=LB)
 INIT_PHI= 1;   % Start simulation with a noisy phi and moments
+MU_P    = 0.0;     % Hermite  hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J    = 0.0;     % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
 % Compute processes distribution
 Ntot = NP_P * NP_KX;
 Nnodes = ceil(Ntot/48);
-- 
GitLab