From 7ee703feb0b9b4db40801643cf5d5bbc236d8b76 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 8 Apr 2022 13:51:33 +0200
Subject: [PATCH] xphij arrays and fort.90 + marconi batch script

---
 .gitignore                 |  4 +-
 batch_script.sh            | 13 +++++++
 fort.90                    | 80 ++++++++++++++++++++++++++++++++++++++
 src/array_mod.F90          |  3 +-
 src/memory.F90             | 11 ++++--
 src/moments_eq_rhs_mod.F90 | 12 +++---
 src/numerics_mod.F90       | 36 +++++++++++++----
 src/restarts_mod.F90       |  4 +-
 wk/analysis_3D.m           |  2 +-
 wk/analysis_header.m       |  2 +-
 10 files changed, 143 insertions(+), 24 deletions(-)
 create mode 100644 batch_script.sh
 create mode 100644 fort.90

diff --git a/.gitignore b/.gitignore
index 963a466a..c03a3145 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,7 +34,7 @@ iCa/
 *.out
 src/srcinfo.h
 src/srcinfo/srcinfo.h
-fort*.90
+*/fort*.90
 local/
-*.sh
+*/*.sh
 Gallery/
diff --git a/batch_script.sh b/batch_script.sh
new file mode 100644
index 00000000..34ffd2e7
--- /dev/null
+++ b/batch_script.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+#SBATCH --job-name=HeLaZ
+#SBATCH --time=24:00:00
+#SBATCH --nodes=1
+#SBATCH --cpus-per-task=1
+#SBATCH --ntasks-per-node=48
+#SBATCH --mem=64GB
+#SBATCH --error=err.txt
+#SBATCH --output=out.txt
+#SBATCH --account=FUA36_TSVVT422
+#SBATCH --partition=skl_fua_dbg
+module load autoload hdf5 fftw
+srun --cpu-bind=cores ./../../../bin/helaz3 2 12 2
diff --git a/fort.90 b/fort.90
new file mode 100644
index 00000000..97a73515
--- /dev/null
+++ b/fort.90
@@ -0,0 +1,80 @@
+&BASIC
+  nrun   = 100000000
+  dt     = 0.01
+  tmax   = 200
+  maxruntime = 356400
+/
+&GRID
+  pmaxe  = 4
+  jmaxe  = 2
+  pmaxi  = 4
+  jmaxi  = 2
+  Nx     = 128
+  Lx     = 100
+  Ny     = 128
+  Ly     = 120
+  Nz     = 16
+  q0     = 1.4
+  shear  = 0
+  eps    = 0.18
+  SG     = .t.
+/
+&OUTPUT_PAR
+  nsave_0d = 200
+  nsave_1d = -1
+  nsave_2d = 200
+  nsave_3d = 200
+  nsave_5d = 2000
+  write_doubleprecision = .f.
+  write_gamma = .t.
+  write_hf    = .t.
+  write_phi   = .t.
+  write_Na00  = .t.
+  write_Napj  = .t.
+  write_Sapj  = .f.
+  write_dens  = .t.
+  write_temp  = .t.
+  job2load    = -1
+/
+&MODEL_PAR
+  ! Collisionality
+  CLOS    = 0
+  NL_CLOS = 0
+  LINEARITY = 'nonlinear'
+  KIN_E   = .f.
+  mu_x    = 0.0
+  mu_y    = 0.0
+  mu_z    = 0.2
+  mu_p    = 0
+  mu_j    = 0
+  nu      = 0.01
+  tau_e   = 1
+  tau_i   = 1
+  sigma_e = 0.023338
+  sigma_i = 1
+  q_e     = -1
+  q_i     = 1
+  K_n     = 2.22
+  K_T     = 6.96
+  K_E     = 0
+  GradB     = 1
+  CurvB     = 1
+  lambdaD = 0
+/
+&COLLISION_PAR
+  collision_model = 'DG'
+  gyrokin_CO      = .true.
+  interspecies    = .true.
+!  mat_file        = '../../../iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+  mat_file        = '../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5'
+/
+&INITIAL_CON
+  INIT_OPT    = 'ppj'
+  ACT_ON_MODES       = 'donothing'
+  init_background  = 0
+  init_noiselvl = 0.00005
+  iseed         = 42
+/
+&TIME_INTEGRATION_PAR
+  numerical_scheme = 'RK4'
+/
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index bafdc82e..841f80d2 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -46,7 +46,8 @@ MODULE array
   REAL(dp), DIMENSION(:),   ALLOCATABLE :: xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1
   REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ynipp1j, ynipm1j,   ynipp1jm1, ynipm1jm1 ! mirror lin coeff for non adiab mom
   REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1            ! mirror lin coeff for adiab mom
-  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_e, xphijp1_e, xphijm1_e
+  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij_i, xphijp1_i, xphijm1_i
   ! Geometrical operators
   ! Curvature
   REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky  ! dimensions: kx, ky, z, odd/even p
diff --git a/src/memory.F90 b/src/memory.F90
index 72d8331f..5c94ee7e 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -86,9 +86,14 @@ SUBROUTINE memory
   CALL allocate_array( dnjs, 1,maxj+1, 1,maxj+1, 1,maxj+1)
 
   ! elect. pot. linear terms
-  CALL allocate_array( xphij,   ips_i,ipe_i, ijs_i,ije_i)
-  CALL allocate_array( xphijp1, ips_i,ipe_i, ijs_i,ije_i)
-  CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
+  IF (KIN_E) THEN
+    CALL allocate_array( xphij_e,   ips_e,ipe_e, ijs_e,ije_e)
+    CALL allocate_array( xphijp1_e, ips_e,ipe_e, ijs_e,ije_e)
+    CALL allocate_array( xphijm1_e, ips_e,ipe_e, ijs_e,ije_e)
+  ENDIF
+  CALL allocate_array( xphij_i,   ips_i,ipe_i, ijs_i,ije_i)
+  CALL allocate_array( xphijp1_i, ips_i,ipe_i, ijs_i,ije_i)
+  CALL allocate_array( xphijm1_i, ips_i,ipe_i, ijs_i,ije_i)
 
   ! Curvature and geometry
   CALL allocate_array( Ckxky,   ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 87380472..0934b533 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -99,9 +99,9 @@ SUBROUTINE moments_eq_rhs_e
           Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
           !! Electrical potential term
           IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
-            Tphi = (xphij  (ip,ij)*kernel_e(ij  ,ikx,iky,iz,eo) &
-                  + xphijp1(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
-                  + xphijm1(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
+            Tphi = (xphij_i  (ip,ij)*kernel_e(ij  ,ikx,iky,iz,eo) &
+                  + xphijp1_i(ip,ij)*kernel_e(ij+1,ikx,iky,iz,eo) &
+                  + xphijm1_i(ip,ij)*kernel_e(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
           ELSE
             Tphi = 0._dp
           ENDIF
@@ -227,9 +227,9 @@ SUBROUTINE moments_eq_rhs_i
 
           !! Electrical potential term
           IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
-            Tphi = (xphij  (ip,ij)*kernel_i(ij  ,ikx,iky,iz,eo) &
-                  + xphijp1(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
-                  + xphijm1(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
+            Tphi = (xphij_i  (ip,ij)*kernel_i(ij  ,ikx,iky,iz,eo) &
+                  + xphijp1_i(ip,ij)*kernel_i(ij+1,ikx,iky,iz,eo) &
+                  + xphijm1_i(ip,ij)*kernel_i(ij-1,ikx,iky,iz,eo))*phi(ikx,iky,iz)
           ELSE
             Tphi = 0._dp
           ENDIF
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index d8e2b178..ba8d0653 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -227,6 +227,27 @@ SUBROUTINE compute_lin_coeff
   ENDDO
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! ES linear coefficients for moment RHS !!!!!!!!!!
+  IF (KIN_E) THEN
+    DO ip = ips_e, ipe_e
+      p_int= parray_e(ip)   ! Hermite degree
+      DO ij = ijs_e, ije_e
+        j_int= jarray_e(ij)   ! REALof Laguerre degree
+        j_dp = REAL(j_int,dp) ! REALof Laguerre degree
+        !! Electrostatic potential pj terms
+        IF (p_int .EQ. 0) THEN ! kronecker p0
+          xphij_e(ip,ij)    =+K_n + 2.*j_dp*K_T
+          xphijp1_e(ip,ij)  =-K_T*(j_dp+1._dp)
+          xphijm1_e(ip,ij)  =-K_T* j_dp
+        ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
+          xphij_e(ip,ij)    =+K_T/SQRT2
+          xphijp1_e(ip,ij)  = 0._dp; xphijm1_e(ip,ij)  = 0._dp;
+        ELSE
+          xphij_e(ip,ij)    = 0._dp; xphijp1_e(ip,ij)  = 0._dp
+          xphijm1_e(ip,ij)  = 0._dp;
+        ENDIF
+      ENDDO
+    ENDDO
+  ENDIF
   DO ip = ips_i, ipe_i
     p_int= parray_i(ip)   ! Hermite degree
     DO ij = ijs_i, ije_i
@@ -234,19 +255,18 @@ SUBROUTINE compute_lin_coeff
       j_dp = REAL(j_int,dp) ! REALof Laguerre degree
       !! Electrostatic potential pj terms
       IF (p_int .EQ. 0) THEN ! kronecker p0
-        xphij(ip,ij)    =+K_n + 2.*j_dp*K_T
-        xphijp1(ip,ij)  =-K_T*(j_dp+1._dp)
-        xphijm1(ip,ij)  =-K_T* j_dp
+        xphij_i(ip,ij)    =+K_n + 2.*j_dp*K_T
+        xphijp1_i(ip,ij)  =-K_T*(j_dp+1._dp)
+        xphijm1_i(ip,ij)  =-K_T* j_dp
       ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
-        xphij(ip,ij)    =+K_T/SQRT2
-        xphijp1(ip,ij)  = 0._dp; xphijm1(ip,ij)  = 0._dp;
+        xphij_i(ip,ij)    =+K_T/SQRT2
+        xphijp1_i(ip,ij)  = 0._dp; xphijm1_i(ip,ij)  = 0._dp;
       ELSE
-        xphij(ip,ij)    = 0._dp; xphijp1(ip,ij)  = 0._dp
-        xphijm1(ip,ij)  = 0._dp;
+        xphij_i(ip,ij)    = 0._dp; xphijp1_i(ip,ij)  = 0._dp
+        xphijm1_i(ip,ij)  = 0._dp;
       ENDIF
     ENDDO
   ENDDO
-
 END SUBROUTINE compute_lin_coeff
 
 !******************************************************************************!
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index cc69a2b1..613cf37f 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -72,10 +72,10 @@ CONTAINS
           ! Read state of system from checkpoint file
           IF (KIN_E) THEN
           WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
-          CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
+          CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
           ENDIF
           WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
-          CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3/))
+          CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikxs:ikxe, ikys:ikye, izs:ize, 1),(/1,3,5/))
 
           CALL closef(fidrst)
 
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index a6f008af..c7a82f47 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -11,7 +11,7 @@ system(['mkdir -p ',MISCDIR]);
 CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
 system(CMD);
 % Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 10; 
+JOBNUMMIN = 00; JOBNUMMAX = 0; 
 data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
 
 
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index a2017024..657b52a8 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -3,6 +3,6 @@
 outfile ='';
 outfile ='linear_shearless_cyclone/dbg';
 % outfile ='debug/test_zpar/';
-% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0_colless/';
+outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0/';
 % outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_0.8/';
 analysis_3D
\ No newline at end of file
-- 
GitLab