From 088f0a5749853b662c3616e6febed284af0a39c3 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Sun, 12 Mar 2023 12:05:26 +0100
Subject: [PATCH] save last changes
---
Makefile | 2 +-
scripts/job_pipeline.sh | 2 +-
src/advance_field_mod.F90 | 4 +-
src/control.F90 | 18 +--
src/diagnose.F90 | 44 +++---
src/geometry_mod.F90 | 170 +++++++++++----------
src/grid_mod.F90 | 159 +++++++++----------
src/moments_eq_rhs_mod.F90 | 42 ++---
src/numerics_mod.F90 | 86 +++++------
src/processing_mod.F90 | 33 ++--
src/solve_EM_fields.F90 | 9 +-
src/tesend.F90 | 2 +-
testcases/smallest_problem/fort.90 | 8 +-
testcases/smallest_problem/fort_00.90 | 8 +-
testcases/smallest_problem/gyacomo_1_debug | 2 +-
wk/header_3D_results.m | 4 +-
16 files changed, 295 insertions(+), 298 deletions(-)
diff --git a/Makefile b/Makefile
index 145be064..34c63b40 100644
--- a/Makefile
+++ b/Makefile
@@ -28,7 +28,7 @@ fast: F90FLAGS = -fast
fast: $(EFST)
# Debug version with all flags
debug: dirs src/srcinfo.h
-debug: F90FLAGS = -C -g -traceback -ftrapuv -warn all -debug all
+debug: F90FLAGS = -g -traceback -ftrapuv -warn all -debug all
# debug: F90FLAGS = -g -traceback -check all -ftrapuv -warn all -debug all
debug: $(EDBG)
# Alpha version, optimized as all but creates another binary
diff --git a/scripts/job_pipeline.sh b/scripts/job_pipeline.sh
index 7cf4f334..342eaa46 100644
--- a/scripts/job_pipeline.sh
+++ b/scripts/job_pipeline.sh
@@ -25,7 +25,7 @@ for i in {1..1}; do
else if (NR == 12) print "srun --cpu-bind=cores ./gyacomo 2 24 1 "ID;
else print $0}' submit_$idm1.cmd > submit_$id.cmd
- # Create new fort file from older one
+ # Create new input file from older one
awk -v "NU=$nu_" -v "TM=$Tm_" -v "J2L=$im1" '{
if (NR == 04) print " tmax = "TM;
else if (NR == 40) print " job2load = "J2L;
diff --git a/src/advance_field_mod.F90 b/src/advance_field_mod.F90
index fd203489..5bed8086 100644
--- a/src/advance_field_mod.F90
+++ b/src/advance_field_mod.F90
@@ -38,7 +38,7 @@ CONTAINS
DO ip =1,local_np
ipi = ip+ngp/2
DO ia =1,local_na
- IF((CLOS .NE. 1) .OR. (parray(ip)+2*jarray(ij) .LE. dmax))&
+ IF((CLOS .NE. 1) .OR. (parray(ipi)+2*jarray(iji) .LE. dmax))&
moments(ia,ipi,iji,iky,ikx,izi,1) = moments(ia,ipi,iji,iky,ikx,izi,1) &
+ dt*b_E(istage)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
END DO
@@ -60,7 +60,7 @@ CONTAINS
DO ip =1,local_np
ipi = ip+ngp/2
DO ia =1,local_na
- IF((CLOS .NE. 1) .OR. (parray(ip)+2*jarray(ij) .LE. dmax))&
+ IF((CLOS .NE. 1) .OR. (parray(ipi)+2*jarray(iji) .LE. dmax))&
moments(ia,ipi,iji,iky,ikx,izi,updatetlevel) = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel) + &
dt*A_E(updatetlevel,istage)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
END DO
diff --git a/src/control.F90 b/src/control.F90
index 74e41c9d..df0bee29 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -60,19 +60,19 @@ SUBROUTINE control
!________________________________________________________________________________
! 2. Main loop
DO
- CALL cpu_time(t0_step) ! Measuring time
-
- CALL tesend
- IF( nlend ) EXIT ! exit do loop
+ CALL cpu_time(t0_step) ! Measuring time
+
+ CALL tesend
+ IF( nlend ) EXIT ! exit do loop
- CALL increase_step
- CALL increase_cstep
- CALL stepon
+ CALL increase_step
+ CALL increase_cstep
- CALL increase_time
+ CALL stepon
+ CALL increase_time
- CALL diagnose(step)
+ CALL diagnose(step)
CALL cpu_time(t1_step);
tc_step = tc_step + (t1_step - t0_step)
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index f1a7d8ee..00ec221c 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -140,22 +140,22 @@ SUBROUTINE diagnose_full(kstep)
CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
! Metric info
CALL creatg(fidres, "/data/metric", "Metric data")
- CALL putarrnd(fidres, "/data/metric/gxx", gxx(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gxy", gxy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gxz", gxz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gyy", gyy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gyz", gyz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gzz", gzz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatR", hatR(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatB", hatB(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdx", dBdx(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdy", dBdy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdz", dBdz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 3/))
- CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,1), (/1, 1, 2, 4/))
+ CALL putarrnd(fidres, "/data/metric/gxx", gxx((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxy", gxy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxz", gxz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyy", gyy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyz", gyz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gzz", gzz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatR", hatR((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatZ", hatZ((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatB", hatB((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdx", dBdx((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdy", dBdy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdz", dBdz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 3/))
+ CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,(1+ngj/2):(local_nj+ngj/2),1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),1), (/1, 1, 2, 4/))
! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN
CALL creatg(fidres, "/data/var0d", "0d profiles")
@@ -346,15 +346,15 @@ SUBROUTINE diagnose_3d
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
! Write current EM fields
- IF (write_phi) CALL write_field3d_kykxz(phi (:,:,1+ngz/2:local_nz+ngz/2), 'phi')
- IF (write_phi.AND.EM) CALL write_field3d_kykxz(psi (:,:,1+ngz/2:local_nz+ngz/2), 'psi')
+ IF (write_phi) CALL write_field3d_kykxz(phi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'phi')
+ IF (write_phi.AND.EM) CALL write_field3d_kykxz(psi (:,:,(1+ngz/2):(local_nz+ngz/2)), 'psi')
IF (write_Na00) THEN
CALL compute_Napjz_spectrum
DO ia=1,local_na
letter_a = name(ia)(1:1)
IF (CONTAINSp0) THEN
! gyrocenter density
- Na00_ = moments(ia,ip0,ij0,:,:,1+ngz/2:local_nz+ngz/2,updatetlevel)
+ Na00_ = moments(ia,ip0,ij0,:,:,(1+ngz/2):(local_nz+ngz/2),updatetlevel)
CALL write_field3d_kykxz(Na00_, 'N'//letter_a//'00')
ENDIF
! <<Napj>x>y spectrum
@@ -470,8 +470,8 @@ SUBROUTINE diagnose_5d
COMPLEX(dp), DIMENSION(total_na,local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
COMPLEX(dp), DIMENSION(total_na,total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name
- field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
- 1:local_nky,1:local_nkx, (1+ngz/2):(local_nz+ngz/2),updatetlevel)
+ field_sub = field(1:total_na,(1+ngp/2):(local_np+ngp/2),((1+ngj/2)):((local_nj+ngj/2)),&
+ 1:local_nky,1:local_nkx, ((1+ngz/2)):((local_nz+ngz/2)),updatetlevel)
field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN
@@ -506,7 +506,7 @@ SUBROUTINE spit_snapshot_check
INQUIRE(file='check_phi', exist=file_exist)
IF( file_exist ) THEN
IF(my_id.EQ. 0) WRITE(*,*) 'Check file found -> gather phi..'
- CALL gather_xyz(phi(:,:,1+Ngz/2:local_nz+Ngz/2), field_to_check,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL gather_xyz(phi(:,:,(1+Ngz/2):(local_nz+Ngz/2)), field_to_check,local_nky,total_nky,total_nkx,local_nz,total_nz)
IF(my_id.EQ. 0) THEN
WRITE(check_filename,'(a16)') 'check_phi.out'
OPEN(fid_check, file=check_filename, form='formatted')
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index cb59aa62..6afb1fed 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -93,7 +93,7 @@ CONTAINS
END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
- USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, Ngz, kxarray, kyarray, set_kparray, Nzgrid, deltaz
+ USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, Ngz, kxarray, kyarray, set_kparray, nzgrid, deltaz
USE basic, ONLY: speak
USE miller, ONLY: set_miller_parameters, get_miller
USE calculus, ONLY: simpson_rule_z
@@ -105,7 +105,7 @@ CONTAINS
INTEGER :: eo,iz,iky,ikx
! Allocate arrays
- CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,Nzgrid)
+ CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,nzgrid)
!
IF( (total_nky .EQ. 1) .AND. (total_nz .EQ. 1)) THEN !1D perp linear run
CALL speak('1D perpendicular geometry')
@@ -136,7 +136,7 @@ CONTAINS
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
CALL set_kparray(gxx,gxy,gyy,hatB)
- DO eo = 1,Nzgrid
+ DO eo = 1,nzgrid
! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = 1,local_nz+Ngz
G1 = gxx(iz,eo)*gyy(iz,eo)-gxy(iz,eo)*gxy(iz,eo)
@@ -160,8 +160,7 @@ CONTAINS
! Gamma_phipar(iz,eo) = G2/G1
ENDDO
ENDDO
-
-
+ !
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
!
@@ -175,14 +174,14 @@ CONTAINS
!
SUBROUTINE eval_salpha_geometry
- USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z
INTEGER :: iz, eo
alpha_MHD = 0._dp
- DO eo = 1,Nzgrid
+ DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
@@ -228,13 +227,13 @@ CONTAINS
!
SUBROUTINE eval_zpinch_geometry
- USE grid, ONLY : local_nz,Ngz,zarray,Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
implicit none
REAL(dp) :: z
INTEGER :: iz, eo
alpha_MHD = 0._dp
- DO eo = 1,Nzgrid
+ DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
@@ -278,12 +277,12 @@ CONTAINS
!--------------------------------------------------------------------------------
! NOT TESTED
subroutine eval_1D_geometry
- USE grid, ONLY : local_nz,Ngz,zarray, Nzgrid
+ USE grid, ONLY : local_nz,Ngz,zarray, nzgrid
! evaluate 1D perp geometry model
implicit none
REAL(dp) :: z
INTEGER :: iz, eo
- DO eo = 1,Nzgrid
+ DO eo = 1,nzgrid
DO iz = 1,local_nz+Ngz
z = zarray(iz,eo)
@@ -310,13 +309,14 @@ CONTAINS
!
SUBROUTINE set_ikx_zBC_map
- USE grid, ONLY: local_nky,Nkx, contains_zmin,contains_zmax, Nexc
+ USE grid, ONLY: local_nky,total_nkx,contains_zmin,contains_zmax, Nexc,&
+ local_nky_offset
USE prec_const, ONLY: imagu, pi
IMPLICIT NONE
! REAL(dp) :: shift
- INTEGER :: ikx,iky
- ALLOCATE(ikx_zBC_L(local_nky,Nkx))
- ALLOCATE(ikx_zBC_R(local_nky,Nkx))
+ INTEGER :: ikx,iky, mn_y
+ ALLOCATE(ikx_zBC_L(local_nky,total_nkx))
+ ALLOCATE(ikx_zBC_R(local_nky,total_nkx))
ALLOCATE(pb_phase_L(local_nky))
ALLOCATE(pb_phase_R(local_nky))
!! No shear case (simple id mapping) or not at the end of the z domain
@@ -326,7 +326,7 @@ CONTAINS
!0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
!(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=free)
DO iky = 1,local_nky
- DO ikx = 1,Nkx
+ DO ikx = 1,total_nkx
ikx_zBC_L(iky,ikx) = ikx ! connect to itself per default
ikx_zBC_R(iky,ikx) = ikx
ENDDO
@@ -340,27 +340,29 @@ CONTAINS
! Modify connection map only at border of z (matters for MPI z-parallelization)
IF(contains_zmin) THEN ! Check if the process is at the start of the fluxtube
DO iky = 1,local_nky
- ! Formula for the shift due to shear after Npol turns
- ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
- DO ikx = 1,Nkx
- ! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc
- ikx_zBC_L(iky,ikx) = ikx-(iky-1)*Nexc
- ! Check if it points out of the kx domain
- ! IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN
- IF( (ikx-(iky-1)*Nexc) .LT. 1 ) THEN ! outside of the frequ domain
- SELECT CASE(parallel_bc)
- CASE ('dirichlet')! connected to 0
- ikx_zBC_L(iky,ikx) = -99
- CASE ('periodic')
- ikx_zBC_L(iky,ikx) = ikx
- CASE ('cyclic')! reroute it by cycling through modes
- ikx_zBC_L(iky,ikx) = MODULO(ikx_zBC_L(iky,ikx)-1,Nkx)+1
- END SELECT
- ENDIF
- ENDDO
- ! phase present in GENE from a shift of the x origin by Lx/2 (useless?)
- ! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_L(iky) = (-1._dp)**(Nexc*(iky-1))*EXP(imagu*REAL(iky-1,dp)*2._dp*pi*shift_y)
+ ! get the real mode number (iky starts at 1 and is shifted from paral)
+ mn_y = iky-1+local_nky_offset
+ ! Formula for the shift due to shear after Npol turns
+ ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
+ DO ikx = 1,total_nkx
+ ! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc
+ ikx_zBC_L(iky,ikx) = ikx-mn_y*Nexc
+ ! Check if it points out of the kx domain
+ ! IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN
+ IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN ! outside of the frequ domain
+ SELECT CASE(parallel_bc)
+ CASE ('dirichlet')! connected to 0
+ ikx_zBC_L(iky,ikx) = -99
+ CASE ('periodic')
+ ikx_zBC_L(iky,ikx) = ikx
+ CASE ('cyclic')! reroute it by cycling through modes
+ ikx_zBC_L(iky,ikx) = MODULO(ikx_zBC_L(iky,ikx)-1,total_nkx)+1
+ END SELECT
+ ENDIF
+ ENDDO
+ ! phase present in GENE from a shift of the x origin by Lx/2 (useless?)
+ ! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
+ pb_phase_L(iky) = (-1._dp)**(Nexc*mn_y)*EXP(imagu*REAL(mn_y,dp)*2._dp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
@@ -368,28 +370,30 @@ CONTAINS
!!!!!!!!!! RIGHT PARALLEL BOUNDARY
IF(contains_zmax) THEN ! Check if the process is at the end of the flux-tube
DO iky = 1,local_nky
- ! Formula for the shift due to shear after Npol
- ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
- DO ikx = 1,Nkx
- ! Usual formula for shifting indices
- ikx_zBC_R(iky,ikx) = ikx+(iky-1)*Nexc
- ! Check if it points out of the kx domain
- ! IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
- IF( (ikx+(iky-1)*Nexc) .GT. Nkx ) THEN ! outside of the frequ domain
- SELECT CASE(parallel_bc)
- CASE ('dirichlet') ! connected to 0
- ikx_zBC_R(iky,ikx) = -99
- CASE ('periodic') ! connected to itself as for shearless
- ikx_zBC_R(iky,ikx) = ikx
- CASE ('cyclic')
- ! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max
- ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,Nkx)+1
- END SELECT
- ENDIF
- ENDDO
- ! phase present in GENE from a shift ofthe x origin by Lx/2 (useless?)
- ! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_R(iky) = (-1._dp)**(Nexc*(iky-1))*EXP(-imagu*REAL(iky-1,dp)*2._dp*pi*shift_y)
+ ! get the real mode number (iky starts at 1 and is shifted from paral)
+ mn_y = iky-1+local_nky_offset
+ ! Formula for the shift due to shear after Npol
+ ! shift = 2._dp*PI*shear*kyarray(iky)*Npol
+ DO ikx = 1,total_nkx
+ ! Usual formula for shifting indices
+ ikx_zBC_R(iky,ikx) = ikx+mn_y*Nexc
+ ! Check if it points out of the kx domain
+ ! IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
+ IF( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain
+ SELECT CASE(parallel_bc)
+ CASE ('dirichlet') ! connected to 0
+ ikx_zBC_R(iky,ikx) = -99
+ CASE ('periodic') ! connected to itself as for shearless
+ ikx_zBC_R(iky,ikx) = ikx
+ CASE ('cyclic')
+ ! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max
+ ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,total_nkx)+1
+ END SELECT
+ ENDIF
+ ENDDO
+ ! phase present in GENE from a shift ofthe x origin by Lx/2 (useless?)
+ ! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
+ pb_phase_R(iky) = (-1._dp)**(Nexc*mn_y)*EXP(-imagu*REAL(mn_y,dp)*2._dp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
@@ -450,31 +454,31 @@ END SUBROUTINE set_ikx_zBC_map
!--------------------------------------------------------------------------------
!
- SUBROUTINE geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,Nzgrid)
- INTEGER, INTENT(IN) :: local_nky,local_nkx,local_nz,Ngz,Nzgrid
+ SUBROUTINE geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,nzgrid)
+ INTEGER, INTENT(IN) :: local_nky,local_nkx,local_nz,Ngz,nzgrid
! Curvature and geometry
- ALLOCATE( Ckxky(local_nky,local_nkx,local_nz+Ngz,Nzgrid))
- ALLOCATE( Jacobian(local_nz+Ngz,Nzgrid))
- ALLOCATE( gxx(local_nz+Ngz,Nzgrid))
- ALLOCATE( gxy(local_nz+Ngz,Nzgrid))
- ALLOCATE( gxz(local_nz+Ngz,Nzgrid))
- ALLOCATE( gyy(local_nz+Ngz,Nzgrid))
- ALLOCATE( gyz(local_nz+Ngz,Nzgrid))
- ALLOCATE( gzz(local_nz+Ngz,Nzgrid))
- ALLOCATE( dBdx(local_nz+Ngz,Nzgrid))
- ALLOCATE( dBdy(local_nz+Ngz,Nzgrid))
- ALLOCATE( dBdz(local_nz+Ngz,Nzgrid))
- ALLOCATE( dlnBdz(local_nz+Ngz,Nzgrid))
- ALLOCATE( hatB(local_nz+Ngz,Nzgrid))
- ! ALLOCATE(Gamma_phipar,(local_nz+Ngz,Nzgrid)) (not implemented)
- ALLOCATE( hatR(local_nz+Ngz,Nzgrid))
- ALLOCATE( hatZ(local_nz+Ngz,Nzgrid))
- ALLOCATE( Rc(local_nz+Ngz,Nzgrid))
- ALLOCATE( phic(local_nz+Ngz,Nzgrid))
- ALLOCATE( Zc(local_nz+Ngz,Nzgrid))
- ALLOCATE( dxdR(local_nz+Ngz,Nzgrid))
- ALLOCATE( dxdZ(local_nz+Ngz,Nzgrid))
- ALLOCATE(gradz_coeff(local_nz+Ngz,Nzgrid))
+ ALLOCATE( Ckxky(local_nky,local_nkx,local_nz+Ngz,nzgrid))
+ ALLOCATE( Jacobian(local_nz+Ngz,nzgrid))
+ ALLOCATE( gxx(local_nz+Ngz,nzgrid))
+ ALLOCATE( gxy(local_nz+Ngz,nzgrid))
+ ALLOCATE( gxz(local_nz+Ngz,nzgrid))
+ ALLOCATE( gyy(local_nz+Ngz,nzgrid))
+ ALLOCATE( gyz(local_nz+Ngz,nzgrid))
+ ALLOCATE( gzz(local_nz+Ngz,nzgrid))
+ ALLOCATE( dBdx(local_nz+Ngz,nzgrid))
+ ALLOCATE( dBdy(local_nz+Ngz,nzgrid))
+ ALLOCATE( dBdz(local_nz+Ngz,nzgrid))
+ ALLOCATE( dlnBdz(local_nz+Ngz,nzgrid))
+ ALLOCATE( hatB(local_nz+Ngz,nzgrid))
+ ! ALLOCATE(Gamma_phipar,(local_nz+Ngz,nzgrid)) (not implemented)
+ ALLOCATE( hatR(local_nz+Ngz,nzgrid))
+ ALLOCATE( hatZ(local_nz+Ngz,nzgrid))
+ ALLOCATE( Rc(local_nz+Ngz,nzgrid))
+ ALLOCATE( phic(local_nz+Ngz,nzgrid))
+ ALLOCATE( Zc(local_nz+Ngz,nzgrid))
+ ALLOCATE( dxdR(local_nz+Ngz,nzgrid))
+ ALLOCATE( dxdZ(local_nz+Ngz,nzgrid))
+ ALLOCATE(gradz_coeff(local_nz+Ngz,nzgrid))
END SUBROUTINE geometry_allocate_mem
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 5b5794fd..9003f696 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -8,10 +8,10 @@ MODULE grid
PRIVATE
! GRID Input
- INTEGER, PUBLIC, PROTECTED :: pmax = 1 ! The maximal Hermite-moment computed
- INTEGER, PUBLIC, PROTECTED :: jmax = 1 ! The maximal Laguerre-moment computed
+ INTEGER, PUBLIC, PROTECTED :: pmax = 1 ! The maximal Hermite-moment computed
+ INTEGER, PUBLIC, PROTECTED :: jmax = 1 ! The maximal Laguerre-moment computed
INTEGER, PUBLIC, PROTECTED :: maxj = 1 ! The maximal Laguerre-moment
- INTEGER, PUBLIC, PROTECTED :: dmax = 1 ! The maximal full GF set of i-moments v^dmax
+ INTEGER, PUBLIC, PROTECTED :: dmax = 1 ! The maximal full GF set of i-moments v^dmax
INTEGER, PUBLIC, PROTECTED :: Nx = 4 ! Number of total internal grid points in x
REAL(dp), PUBLIC, PROTECTED :: Lx = 120_dp ! horizontal length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nexc = 1 ! factor to increase Lx when shear>0 (Lx = Nexc/kymin/shear)
@@ -19,8 +19,8 @@ MODULE grid
REAL(dp), PUBLIC, PROTECTED :: Ly = 120_dp ! vertical length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nz = 4 ! Number of total perpendicular planes
INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes
- INTEGER, PUBLIC, PROTECTED :: Nkx = 4 ! Number of total internal grid points in kx
- INTEGER, PUBLIC, PROTECTED :: Nky = 4 ! Number of total internal grid points in ky
+ INTEGER, PUBLIC, PROTECTED :: Nkx ! Number of total internal grid points in kx
+ INTEGER, PUBLIC, PROTECTED :: Nky ! Number of total internal grid points in ky
REAL(dp), PUBLIC, PROTECTED :: kpar = 0_dp ! parallel wave vector component
! Grid arrays
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: parray, parray_full
@@ -82,8 +82,8 @@ MODULE grid
REAL(dp), PUBLIC, PROTECTED :: diff_kx_coeff, diff_ky_coeff, diff_dz_coeff
LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag
! Array to know the distribution of data among all processes (for MPI comm)
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: counts_nkx, counts_nky, counts_nz
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: displs_nkx, displs_nky, displs_nz
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: counts_total_nkx, counts_nky, counts_nz
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: displs_total_nkx, displs_nky, displs_nz
! Kperp array depends on kx, ky, z (geometry), eo (even or odd zgrid)
LOGICAL, PUBLIC, PROTECTED :: contains_kx0 = .false. ! flag if the proc contains kx=0 index
LOGICAL, PUBLIC, PROTECTED :: contains_ky0 = .false. ! flag if the proc contains ky=0 index
@@ -164,11 +164,15 @@ CONTAINS
CALL set_kxgrid(shear,Npol,LINEARITY,N_HD)
CALL set_zgrid (Npol)
- print*, 'p:',parray
- print*, 'j:',jarray
- print*, 'ky:',kyarray
- print*, 'kx:',kxarray
- print*, 'z:',zarray
+ ! print*, 'p:',parray
+ ! print*, 'j:',jarray
+ ! print*, 'ky:',kyarray
+ ! print*, 'kx:',kxarray
+ ! print*, 'z:',zarray
+ ! print*, parray(ip0)
+ ! print*, jarray(ij0)
+ ! print*, kyarray(iky0)
+ ! print*, kxarray(ikx0)
END SUBROUTINE set_grids
SUBROUTINE init_1Dgrid_distr
@@ -303,7 +307,7 @@ CONTAINS
IMPLICIT NONE
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
- INTEGER :: iky, ikyo
+ INTEGER :: iky
Nky = Ny/2+1 ! Defined only on positive kx since fields are real
! Grid spacings
IF (Ny .EQ. 1) THEN
@@ -329,36 +333,35 @@ CONTAINS
local_kymax = 0._dp
! Creating a grid ordered as dk*(0 1 2 3)
! We loop over the natural iky numbers (|1 2 3||4 5 6||... Nky|)
- DO iky = ikys,ikye
+ DO iky = 1,local_nky
! We shift the natural iky index by the offset to obtain the mpi dependent
- ! indexation (|1 2 3||1 2 3|... local_Nky|)
- ikyo = iky - local_nky_offset
+ ! indexation (|1 2 3||1 2 3|... local_nky|)
IF(Ny .EQ. 1) THEN
kyarray(iky) = deltaky
kyarray_full(iky) = deltaky
SINGLE_KY = .TRUE.
ELSE
- kyarray(ikyo) = REAL(iky,dp) * deltaky
+ kyarray(iky) = kyarray_full(iky-local_nky_offset)
ENDIF
! Finding kx=0
- IF (kyarray(ikyo) .EQ. 0) THEN
- iky0 = ikyo
+ IF (kyarray(iky) .EQ. 0) THEN
+ iky0 = iky
contains_ky0 = .true.
ENDIF
! Finding local kxmax value
- IF (ABS(kyarray(ikyo)) .GT. local_kymax) THEN
- local_kymax = ABS(kyarray(ikyo))
+ IF (ABS(kyarray(iky)) .GT. local_kymax) THEN
+ local_kymax = ABS(kyarray(iky))
ENDIF
! Finding kxmax idx
- IF (kyarray(ikyo) .EQ. ky_max) THEN
- iky_max = ikyo
+ IF (kyarray(iky) .EQ. ky_max) THEN
+ iky_max = iky
contains_kymax = .true.
ENDIF
END DO
! Orszag 2/3 filter
two_third_kymax = 2._dp/3._dp*deltaky*(Nky-1)
- ALLOCATE(AA_y(local_Nky))
- DO iky = 1,local_Nky
+ ALLOCATE(AA_y(local_nky))
+ DO iky = 1,local_nky
IF ( (kyarray(iky) .LT. two_third_kymax) .OR. (LINEARITY .EQ. 'linear')) THEN
AA_y(iky) = 1._dp;
ELSE
@@ -380,7 +383,7 @@ CONTAINS
INTEGER, INTENT(IN) :: Npol
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
- INTEGER :: ikx, ikxo
+ INTEGER :: ikx
REAL(dp):: Lx_adapted
IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
@@ -394,17 +397,17 @@ CONTAINS
! x length is adapted
Lx = Lx_adapted*Nexc
ENDIF
- Nkx = Nx;
+ Nkx = Nx
total_nkx = Nx
! Local data
! Start and END indices of grid
ikxs = 1
- ikxe = Nkx
+ ikxe = total_nkx
local_nkx_ptr = ikxe - ikxs + 1
local_nkx = ikxe - ikxs + 1
local_nky_offset = ikxs - 1
ALLOCATE(kxarray(local_nkx))
- ALLOCATE(kxarray_full(1:total_nkx))
+ ALLOCATE(kxarray_full(total_nkx))
IF (Nx .EQ. 1) THEN
deltakx = 1._dp
kxarray(1) = 0._dp
@@ -417,66 +420,31 @@ CONTAINS
local_kxmax = 0._dp
ELSE ! Build apprpopriate grid
deltakx = 2._dp*PI/Lx
- IF(MODULO(Nkx,2) .EQ. 0) THEN ! Even number of Nkx (-2 -1 0 1 2 3)
- kx_max = (Nkx/2)*deltakx
+ IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
+ kx_max = (total_nkx/2)*deltakx
kx_min = -kx_max+deltakx
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
- local_kxmax = 0._dp
- DO ikx = ikxs,ikxe
- ikxo = ikx - local_nkx_offset
- kxarray(ikxo) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
- if (ikx .EQ. Nx/2+1) kxarray(ikxo) = -kxarray(ikxo)
- ! Finding kx=0
- IF (kxarray(ikxo) .EQ. 0) THEN
- ikx0 = ikxo
- contains_kx0 = .true.
- ENDIF
- ! Finding local kxmax
- IF (ABS(kxarray(ikxo)) .GT. local_kxmax) THEN
- local_kxmax = ABS(kxarray(ikxo))
- ENDIF
- ! Finding kxmax
- IF (kxarray(ikxo) .EQ. kx_max) ikx_max = ikxo
- END DO
- ! Build the full grids on process 0 to diagnose it without comm
- ! kx
- DO ikx = 1,Nkx
- kxarray_full(ikx) = deltakx*(MODULO(ikx-1,Nkx/2)-Nkx/2*FLOOR(2.*real(ikx-1,dp)/real(Nkx,dp)))
- IF (ikx .EQ. Nx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
+ DO ikx = 1,total_nkx
+ kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),dp)
+ IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
- ELSE ! Odd number of kx (-2 -1 0 1 2)
- kx_max = (Nkx-1)/2*deltakx
- kx_min = -kx_max
- ! Creating a grid ordered as dk*(0 1 2 -2 -1)
+ ! Set local grid (not parallelized so same as full one)
local_kxmax = 0._dp
- DO ikx = ikxs,ikxe
- ikxo = ikx - local_nkx_offset
- IF(ikx .LE. (Nkx-1)/2+1) THEN
- kxarray(ikxo) = deltakx*(ikx-1)
- ELSE
- kxarray(ikxo) = deltakx*(ikx-Nkx-1)
- ENDIF
+ DO ikx = 1,local_nkx
+ kxarray(ikx) = kxarray_full(ikx-local_nkx_offset)
! Finding kx=0
- IF (kxarray(ikxo) .EQ. 0) THEN
- ikx0 = ikxo
+ IF (kxarray(ikx) .EQ. 0) THEN
+ ikx0 = ikx
contains_kx0 = .true.
ENDIF
! Finding local kxmax
- IF (ABS(kxarray(ikxo)) .GT. local_kxmax) THEN
- local_kxmax = ABS(kxarray(ikxo))
- ENDIF
- ! Finding kxmax
- IF (kxarray(ikxo) .EQ. kx_max) ikx_max = ikxo
- END DO
- ! Build the full grids on process 0 to diagnose it without comm
- ! kx
- DO ikx = 1,Nkx
- IF(ikx .LE. (Nkx-1)/2+1) THEN
- kxarray_full(ikx) = deltakx*(ikx-1)
- ELSE
- kxarray_full(ikx) = deltakx*(ikx-Nkx-1)
+ IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
+ local_kxmax = ABS(kxarray(ikx))
+ ikx_max = ikx
ENDIF
END DO
+ ELSE ! Odd number of kx (-2 -1 0 1 2)
+ kx_max = (total_nkx-1)/2*deltakx
ENDIF
ENDIF
! Orszag 2/3 filter
@@ -504,7 +472,7 @@ CONTAINS
USE parallel, ONLY: num_procs_z, rank_z
IMPLICIT NONE
REAL(dp):: grid_shift, Lz, zmax, zmin
- INTEGER :: istart, iend, in, Npol, iz, ig, eo
+ INTEGER :: istart, iend, in, Npol, iz, ig, eo, iglob
total_nz = Nz
! Length of the flux tube (in ballooning angle)
Lz = 2_dp*pi*Npol
@@ -572,11 +540,24 @@ CONTAINS
! Find local extrema
local_zmax(eo) = zarray(local_nz+ngz/2,eo)
local_zmin(eo) = zarray(1+ngz/2,eo)
- print*, zarray
! Fill the ghosts
- DO ig = 1,ngz/2
- zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),dp)*deltaz
- zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,dp)*deltaz
+ ! Continue angles
+ ! DO ig = 1,ngz/2
+ ! zarray(ig,eo) = local_zmin(eo)-REAL(ngz/2-(ig-1),dp)*deltaz
+ ! zarray(local_nz+ngz/2+ig,eo) = local_zmax(eo)+REAL(ig,dp)*deltaz
+ ! ENDDO
+ ! Periodic z \in (-pi pi-dz)
+ DO ig = 1,ngz/2 ! first ghost cells
+ iglob = ig+local_nz_offset-ngz/2
+ IF (iglob .LE. 0) &
+ iglob = iglob + total_nz
+ zarray(ig,eo) = zarray_full(iglob)
+ ENDDO
+ DO ig = local_nz+ngz/2,local_nz+ngz ! last ghost cells
+ iglob = ig+local_nz_offset-ngz/2
+ IF (iglob .GT. total_nz) &
+ iglob = iglob - total_nz
+ zarray(ig,eo) = zarray_full(iglob)
ENDDO
! Set up the flags to know if the process contains the tip and/or the tail
! of the z domain (important for z-boundary condition)
@@ -591,12 +572,12 @@ CONTAINS
END SUBROUTINE set_zgrid
SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
- REAL(dp), DIMENSION(:,:), INTENT(IN) :: gxx,gxy,gyy,hatB
+ REAL(dp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,hatB
INTEGER :: eo,iz,iky,ikx
REAL(dp) :: kx, ky
- CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+Ngz, 1,2)
- DO eo = 1,Nzgrid
- DO iz = 1,local_nz+Ngz
+ CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
+ DO eo = 1,nzgrid
+ DO iz = 1,local_nz+ngz
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx = 1,local_nkx
@@ -628,7 +609,7 @@ CONTAINS
CALL attach(fid, TRIM(str), "Ny", Ny)
CALL attach(fid, TRIM(str), "Ly", Ly)
CALL attach(fid, TRIM(str), "Nz", Nz)
- CALL attach(fid, TRIM(str), "Nkx", Nkx)
+ CALL attach(fid, TRIM(str), "total_nkx", total_nkx)
CALL attach(fid, TRIM(str), "Nky", Nky)
CALL attach(fid, TRIM(str), "SG", SG)
END SUBROUTINE grid_outputinputs
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index db3a6c85..7a127f53 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -43,15 +43,14 @@ SUBROUTINE compute_moments_eq_rhs
ky = kyarray(iky) ! binormal wavevector
i_ky = imagu * ky ! binormal derivative
! Kinetic loops
- j:DO ij = 1, local_nj ! This loop is from 1 to jmaxi+1
+ j:DO ij = 1,local_nj ! This loop is from 1 to jmaxi+1
iji = ij+ngj/2
j_int = jarray(iji)
- p:DO ip = 1, local_np ! Hermite loop
+ p:DO ip = 1,local_np ! Hermite loop
ipi = ip+ngp/2
p_int = parray(ipi) ! Hermite degree
eo = min(nzgrid,MODULO(p_int,2)+1) ! Indicates if we are on odd or even z grid
kperp2= kparray(iky,ikx,izi,eo)**2
- RHS = 0._dp
! Species loop
a:DO ia = 1,local_na
Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
@@ -87,7 +86,7 @@ SUBROUTINE compute_moments_eq_rhs
Fmir = dlnBdz(iz,eo)*(Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 +&
Unapm1j + Unapm1jp1 + Unapm1jm1)
! Parallel magnetic term (Landau damping and the mirror force)
- Mpara = gradz_coeff(iz,eo)*(Ldamp + Fmir)
+ Mpara = gradz_coeff(iz,eo)*(Ldamp) !+ Fmir)
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
Dphi =i_ky*( xphij (ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
@@ -111,20 +110,20 @@ SUBROUTINE compute_moments_eq_rhs
! Perpendicular magnetic term
- Mperp &
! Parallel magnetic term
- - Mpara &
- ! Drives (density + temperature gradients)
- - (Dphi + Dpsi) &
- ! Collision term
- + Capj(ia,ip,ij,iky,ikx,iz) &
- ! Perpendicular pressure effects (electromagnetic term) (TO CHECK)
- - i_ky*beta*dpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
- ! Parallel drive term (should be negligible, to test)
- ! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
- ! Numerical perpendicular hyperdiffusion
- -mu_x*diff_kx_coeff*kx**N_HD*Napj &
- -mu_y*diff_ky_coeff*ky**N_HD*Napj &
- ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
- -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,iz)
+ - Mpara !&
+ ! ! Drives (density + temperature gradients)
+ ! - (Dphi + Dpsi) &
+ ! ! Collision term
+ ! + Capj(ia,ip,ij,iky,ikx,iz) &
+ ! ! Perpendicular pressure effects (electromagnetic term) (TO CHECK)
+ ! - i_ky*beta*dpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
+ ! ! Parallel drive term (should be negligible, to test)
+ ! ! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
+ ! ! Numerical perpendicular hyperdiffusion
+ ! -mu_x*diff_kx_coeff*kx**N_HD*Napj &
+ ! -mu_y*diff_ky_coeff*ky**N_HD*Napj &
+ ! ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
+ ! -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,iz)
!! Velocity space dissipation (should be implemented somewhere else)
SELECT CASE(HYP_V)
CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
@@ -134,7 +133,7 @@ SUBROUTINE compute_moments_eq_rhs
RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
CASE('dvpar4')
! fourth order numerical diffusion in vpar
- IF(ip-4 .GT. 0) &
+ IF(p_int .GE. 4) &
! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
! (not used often so not parallelized)
RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
@@ -157,6 +156,11 @@ SUBROUTINE compute_moments_eq_rhs
! Execution time end
CALL cpu_time(t1_rhs)
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
+ ! print*, moments(1,3,2,2,2,3,updatetlevel)
+ ! print*, moments_rhs(1,3,2,2,2,3,updatetlevel)
+ print*, SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
+ stop
+
END SUBROUTINE compute_moments_eq_rhs
! SUBROUTINE add_Maxwellian_background_terms
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index ff0949bb..98e08e18 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -72,7 +72,7 @@ END SUBROUTINE build_dv4Hp_table
SUBROUTINE evaluate_kernels
USE basic
USE array, ONLY : kernel!, HF_phi_correction_operator
- USE grid, ONLY : local_Na, local_Nj,Ngj, local_nkx, local_nky, local_nz, Ngz, jarray, kparray,&
+ USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kparray,&
nzgrid
USE species, ONLY : sigma2_tau_o2
USE prec_const, ONLY: dp
@@ -80,24 +80,22 @@ SUBROUTINE evaluate_kernels
INTEGER :: j_int, ia, eo, ikx, iky, iz, ij
REAL(dp) :: j_dp, y_, factj
-DO ia = 1,local_Na
+DO ia = 1,local_na
DO eo = 1,nzgrid
DO ikx = 1,local_nkx
DO iky = 1,local_nky
- DO iz = 1,local_nz + Ngz
- DO ij = 1, local_nj + Ngj
+ DO iz = 1,local_nz + ngz
+ DO ij = 1,local_nj + ngj
j_int = jarray(ij)
j_dp = REAL(j_int,dp)
y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
- IF(j_int .LT. 0) THEN
+ IF(j_int .LT. 0) THEN !ghosts values
kernel(ia,ij,iky,ikx,iz,eo) = 0._dp
ELSE
- factj = GAMMA(j_dp+1._dp)
+ factj = REAL(GAMMA(j_dp+1._dp),dp)
kernel(ia,ij,iky,ikx,iz,eo) = y_**j_int*EXP(-y_)/factj
ENDIF
ENDDO
- ! IF (ijs .EQ. 1) & ! if ijs is 1, the ghost kernel has negative index
- ! kernel(ia,ijgs,iky,ikx,iz,eo) = 0._dp
ENDDO
ENDDO
ENDDO
@@ -107,7 +105,7 @@ DO ia = 1,local_Na
! 2._dp * Kernel(ia,1,:,:,:,1) &
! -1._dp * Kernel(ia,2,:,:,:,1)
!
- ! DO ij = 1, local_Nj
+ ! DO ij = 1,local_Nj
! j_int = jarray(ij)
! j_dp = REAL(j_int,dp)
! HF_phi_correction_operator(:,:,:) = HF_phi_correction_operator(:,:,:) &
@@ -117,7 +115,6 @@ DO ia = 1,local_Na
! - j_dp * Kernel(ia,ij-1,:,:,:,1))
! ENDDO
ENDDO
-
END SUBROUTINE evaluate_kernels
!******************************************************************************!
@@ -134,7 +131,7 @@ END SUBROUTINE evaluate_EM_op
SUBROUTINE evaluate_poisson_op
USE basic
USE array, ONLY : kernel, inv_poisson_op, inv_pol_ion
- USE grid, ONLY : local_Na, local_nkx, local_nky, local_nz,&
+ USE grid, ONLY : local_na, local_nkx, local_nky, local_nz,&
kxarray, kyarray, local_nj, ngj, ngz, ieven
USE species, ONLY : q2_tau
USE model, ONLY : ADIAB_E, tau_e
@@ -153,14 +150,14 @@ SUBROUTINE evaluate_poisson_op
ELSE
operator = 0._dp
DO ia = 1,local_na ! sum over species
- pol_tot = 0._dp ! total polarisation term
- pol_ion = 0._dp ! sum of ion polarisation term
! loop over n only up to the max polynomial degree
+ pol_tot = 0._dp ! total polarisation term
+ pol_ion = 0._dp ! sum of ion polarisation term
DO in=1,local_nj
- pol_tot = pol_tot + q2_tau(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 ! ... sum recursively ...
- pol_ion = pol_ion + q2_tau(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 !
+ pol_tot = pol_tot + kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 ! ... sum recursively ...
+ pol_ion = pol_ion + kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)**2 !
END DO
- operator = operator + q2_tau(ia) - pol_tot
+ operator = operator + q2_tau(ia)*(1._dp - pol_tot)
ENDDO
operator_ion = operator
IF(ADIAB_E) THEN ! Adiabatic model
@@ -172,7 +169,6 @@ SUBROUTINE evaluate_poisson_op
END DO zloop
END DO kyloop
END DO kxloop
-
END SUBROUTINE evaluate_poisson_op
!******************************************************************************!
@@ -182,7 +178,7 @@ END SUBROUTINE evaluate_poisson_op
SUBROUTINE evaluate_ampere_op
USE prec_const, ONLY : dp
USE array, ONLY : kernel, inv_ampere_op
- USE grid, ONLY : local_Na, local_nkx, local_nky, local_nz, &
+ USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, &
jmax, kparray, kxarray, kyarray, SOLVE_AMPERE
USE model, ONLY : beta
USE species, ONLY : q, sigma
@@ -226,7 +222,7 @@ SUBROUTINE compute_lin_coeff
xnapp1j, xnapm1j, xnapp2j, xnapm2j,&
xphij, xphijp1, xphijm1,&
xpsij, xpsijp1, xpsijm1
- USE species, ONLY: k_T, k_N, tau, q, sqrtTau_q, tau_q
+ USE species, ONLY: k_T, k_N, tau, q, sqrt_tau_o_sigma
USE model, ONLY: k_cB, k_gB
USE prec_const, ONLY: dp, SQRT2, SQRT3
USE grid, ONLY: parray, jarray, local_na, local_np, local_nj, ngj, ngp
@@ -235,48 +231,48 @@ SUBROUTINE compute_lin_coeff
!! linear coefficients for moment RHS !!!!!!!!!!
DO ia = 1,local_na
- DO ip = 1, local_np
+ DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- DO ij = 1, local_nj
+ DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! All Napj terms
xnapj(ia,ip,ij) = tau(ia)/q(ia)*(k_cB*(2._dp*p_dp + 1._dp) &
- +k_gB*(2._dp*j_dp + 1._dp))
+ +k_gB*(2._dp*j_dp + 1._dp))
! Mirror force terms
- ynapp1j (ia,ip,ij) = -sqrtTau_q(ia) * (j_dp+1._dp)*SQRT(p_dp+1._dp)
- ynapm1j (ia,ip,ij) = -sqrtTau_q(ia) * (j_dp+1._dp)*SQRT(p_dp)
- ynapp1jm1(ia,ip,ij) = +sqrtTau_q(ia) * j_dp*SQRT(p_dp+1._dp)
- ynapm1jm1(ia,ip,ij) = +sqrtTau_q(ia) * j_dp*SQRT(p_dp)
+ ynapp1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp+1._dp)
+ ynapm1j (ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp)
+ ynapp1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp+1._dp)
+ ynapm1jm1(ia,ip,ij) = +sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp)
! Trapping terms
- zNapm1j (ia,ip,ij) = +sqrtTau_q(ia) *(2._dp*j_dp+1._dp)*SQRT(p_dp)
- zNapm1jp1(ia,ip,ij) = -sqrtTau_q(ia) * (j_dp+1._dp)*SQRT(p_dp)
- zNapm1jm1(ia,ip,ij) = -sqrtTau_q(ia) * j_dp*SQRT(p_dp)
+ zNapm1j (ia,ip,ij) = +sqrt_tau_o_sigma(ia) *(2._dp*j_dp+1._dp)*SQRT(p_dp)
+ zNapm1jp1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * (j_dp+1._dp)*SQRT(p_dp)
+ zNapm1jm1(ia,ip,ij) = -sqrt_tau_o_sigma(ia) * j_dp*SQRT(p_dp)
ENDDO
ENDDO
- DO ip = 1, local_np
+ DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
p_dp = REAL(p_int,dp) ! REAL of Hermite degree
! Landau damping coefficients (ddz napj term)
- xnapp1j(ia,ip) = sqrtTau_q(ia) * SQRT(p_dp+1._dp)
- xnapm1j(ia,ip) = sqrtTau_q(ia) * SQRT(p_dp)
+ xnapp1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_dp+1._dp)
+ xnapm1j(ia,ip) = sqrt_tau_o_sigma(ia) * SQRT(p_dp)
! Magnetic curvature term
- xnapp2j(ia,ip) = tau_q(ia) * k_cB * SQRT((p_dp+1._dp)*(p_dp + 2._dp))
- xnapm2j(ia,ip) = tau_q(ia) * k_cB * SQRT( p_dp *(p_dp - 1._dp))
+ xnapp2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT((p_dp+1._dp)*(p_dp + 2._dp))
+ xnapm2j(ia,ip) = tau(ia)/q(ia) * k_cB * SQRT( p_dp *(p_dp - 1._dp))
ENDDO
- DO ij = 1, local_nj
+ DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! Laguerre degree
j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
! Magnetic gradient term
- xnapjp1(ia,ij) = -tau_q(ia) * k_gB * (j_dp + 1._dp)
- xnapjm1(ia,ij) = -tau_q(ia) * k_gB * j_dp
+ xnapjp1(ia,ij) = -tau(ia)/q(ia) * k_gB * (j_dp + 1._dp)
+ xnapjm1(ia,ij) = -tau(ia)/q(ia) * k_gB * j_dp
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ES linear coefficients for moment RHS !!!!!!!!!!
- DO ip = 1, local_np
+ DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
- DO ij = 1, local_nj
+ DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
@@ -295,17 +291,17 @@ SUBROUTINE compute_lin_coeff
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Electromagnatic linear coefficients for moment RHS !!!!!!!!!!
- DO ip = 1, local_np
+ DO ip = 1,local_np
p_int= parray(ip+ngp/2) ! Hermite degree
- DO ij = 1, local_nj
+ DO ij = 1,local_nj
j_int= jarray(ij+ngj/2) ! REALof Laguerre degree
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij (ia,ip,ij) = +(k_N(ia) + (2._dp*j_dp+1._dp)*k_T(ia))* sqrtTau_q(ia)
- xpsijp1(ia,ip,ij) = - k_T(ia)*(j_dp+1._dp) * sqrtTau_q(ia)
- xpsijm1(ia,ip,ij) = - k_T(ia)* j_dp * sqrtTau_q(ia)
+ xpsij (ia,ip,ij) = +(k_N(ia) + (2._dp*j_dp+1._dp)*k_T(ia))* sqrt_tau_o_sigma(ia)
+ xpsijp1(ia,ip,ij) = - k_T(ia)*(j_dp+1._dp) * sqrt_tau_o_sigma(ia)
+ xpsijm1(ia,ip,ij) = - k_T(ia)* j_dp * sqrt_tau_o_sigma(ia)
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij (ia,ip,ij) = + k_T(ia)*SQRT3/SQRT2 * sqrtTau_q(ia)
+ xpsij (ia,ip,ij) = + k_T(ia)*SQRT3/SQRT2 * sqrt_tau_o_sigma(ia)
xpsijp1(ia,ip,ij) = 0._dp; xpsijm1(ia,ip,ij) = 0._dp;
ELSE
xpsij (ia,ip,ij) = 0._dp; xpsijp1(ia,ip,ij) = 0._dp
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 0f88cc6b..5872be31 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -1,7 +1,7 @@
MODULE processing
USE prec_const, ONLY: dp, imagu, SQRT2, SQRT3
USE grid, ONLY: &
- local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp, &
+ local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp,nzgrid, &
parray,pmax,ip0, iodd, ieven,&
CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
jarray,jmax,ij0, dmax,&
@@ -286,32 +286,43 @@ CONTAINS
DO ij = 1,local_nj+ngj
DO ip = 1,local_np+ngp
DO ia = 1,local_na
- p_int = parray(ip)
- eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
+ IF(nzgrid .GT. 1) THEN
+ p_int = parray(ip+ngp/2)
+ eo = MODULO(p_int,2)+1 ! Indicates if we are on even or odd z grid
+ ELSE
+ eo = 0
+ ENDIF
! Compute z first derivative
f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
- ddz_napj(ia,ip,ij,iky,ikx,:) = f_out
+ ! get output
+ DO iz = 1,local_nz
+ ddz_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
+ ENDDO
! Parallel numerical diffusion
IF (HDz_h) THEN
f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
- ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
ELSE
f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
- CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
- ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
ENDIF
+ CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
+ ! get output
+ DO iz = 1,local_nz
+ ddzND_Napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
+ ENDDO
! Compute even odd grids interpolation
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,1:local_nz+ngz)
+ f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
CALL interp_z(eo,local_nz,ngz,f_in,f_out)
- interp_napj(ia,ip,ij,iky,ikx,1:local_nz) = f_out
+ DO iz = 1,local_nz
+ interp_napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
+ ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
-
+ print*, ddz_napj(1,3,2,:,:,:)
+ stop
! Phi parallel gradient (not implemented fully, should be negligible)
! DO ikx = 1,local_nkx
! DO iky = 1,local_nky
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 0b78a8dd..01717c89 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -36,15 +36,14 @@ CONTAINS
CALL cpu_time(t0_poisson)
!! Poisson can be solved only for process containng p=0
IF ( SOLVE_POISSON ) THEN
- x:DO ikx = 1,local_nky
- y:DO iky = 1,local_nkx
+ x:DO ikx = 1,local_nkx
+ y:DO iky = 1,local_nky
!!!!!!!!!!!!!!! Compute particle charge density q_a n_a for each evolved species
DO iz = 1,local_nz
rho(iz) = 0._dp
- DO in=1,total_nj
+ DO in = 1,total_nj
DO ia = 1,local_na
- rho(iz) = rho(iz) &
- +q(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)*moments(ia,ip0,in+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
+ rho(iz) = rho(iz) + q(ia)*kernel(ia,in+ngj/2,iky,ikx,iz+ngz/2,ieven)*moments(ia,ip0,in+ngj/2,iky,ikx,iz+ngz/2,updatetlevel)
END DO
END DO
END DO
diff --git a/src/tesend.F90 b/src/tesend.F90
index 6d53392c..14d67344 100644
--- a/src/tesend.F90
+++ b/src/tesend.F90
@@ -33,7 +33,7 @@ SUBROUTINE tesend
!________________________________________________________________________________
! 3. Test on TMAX
- nlend = time .GT. tmax
+ nlend = time .GE. tmax
IF ( nlend ) THEN
CALL speak('TMAX reached')
RETURN
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index fdb25d1d..8fb8b7da 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -19,7 +19,7 @@
&GEOMETRY
geom = 's-alpha'
q0 = 1.4
- shear = 0.8
+ shear = 0.0
eps = 0.18
kappa = 1.0
s_kappa= 0.0
@@ -77,7 +77,7 @@
name_ = 'electrons'
tau_ = 1.0
sigma_= 0.023338
- q_ = 1.0
+ q_ = -1.0
k_N_ = 0.0!2.22
k_T_ = 0.0!6.96
/
@@ -91,8 +91,8 @@
&INITIAL_CON
INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing'
- init_background = 0.0
- init_noiselvl = 1.0
+ init_background = 1.0
+ init_noiselvl = 0.0
iseed = 42
/
&TIME_INTEGRATION_PAR
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index d37fc3f1..f14e2591 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -19,7 +19,7 @@
&GEOMETRY
geom = 's-alpha'
q0 = 1.4
- shear = 0.8
+ shear = 0.0
eps = 0.18
parallel_bc = 'dirichlet'
/
@@ -44,7 +44,7 @@
CLOS = 0
NL_CLOS = -1
LINEARITY = 'linear'
- KIN_E = .f.
+ KIN_E = .t.
mu_x = 0.0
mu_y = 0.0
N_HD = 4
@@ -73,8 +73,8 @@
&INITIAL_CON
INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing'
- init_background = 0.0
- init_noiselvl = 1.0
+ init_background = 1.0
+ init_noiselvl = 0.0
iseed = 42
/
&TIME_INTEGRATION_PAR
diff --git a/testcases/smallest_problem/gyacomo_1_debug b/testcases/smallest_problem/gyacomo_1_debug
index f5e7015f..f0c8bbfd 120000
--- a/testcases/smallest_problem/gyacomo_1_debug
+++ b/testcases/smallest_problem/gyacomo_1_debug
@@ -1 +1 @@
-../../bin/gyacomo_1_debug
\ No newline at end of file
+../../../gyacomo_1/bin/gyacomo_debug
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 25115ab7..70619eee 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -47,7 +47,9 @@ PARTITION = '/misc/gyacomo_outputs/';
%% kT=5.3 results
% resdir = 'paper_2_nonlinear/kT_5.3/5x3x128x64x64';
-% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_5.3/5x3x128x64x24';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_sp';
+% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x64_dp';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_MUxy_0';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_NL_-1';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_nuDG_0.01';
--
GitLab