diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 83a972fd8e3d7057fc1bcb965193de5883f71e0d..dd60f23bcf700aaeab9e35fa71b2019b2e285157 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -20,12 +20,21 @@ implicit none
REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
REAL(dp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
! parameters for Miller geometry
- REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation
+ REAL(dp), PUBLIC, PROTECTED :: kappa = 1._dp ! elongation (1 for circular)
REAL(dp), PUBLIC, PROTECTED :: s_kappa = 0._dp ! r normalized derivative skappa = r/kappa dkappa/dr
REAL(dp), PUBLIC, PROTECTED :: delta = 0._dp ! triangularity
REAL(dp), PUBLIC, PROTECTED :: s_delta = 0._dp ! '' sdelta = r/sqrt(1-delta2) ddelta/dr
REAL(dp), PUBLIC, PROTECTED :: zeta = 0._dp ! squareness
REAL(dp), PUBLIC, PROTECTED :: s_zeta = 0._dp ! '' szeta = r dzeta/dr
+ ! to apply shift in the parallel z-BC if shearless
+ REAL(dp), PUBLIC, PROTECTED :: shift_y = 0._dp ! for Arno
+ ! Chooses the type of parallel BC we use for the unconnected kx modes (active for non-zero shear only)
+ ! 'periodic' : Connect a disconnected kx to a mode on the other cadran
+ ! 'dirichlet' : Connect a disconnected kx to 0
+ ! 'disconnected' : Connect all kx to 0
+ ! 'shearless' : Connect all kx to itself
+ CHARACTER(len=256), &
+ PUBLIC, PROTECTED :: parallel_bc
! GENE unused additional parameters for miller_mod
REAL(dp), PUBLIC, PROTECTED :: edge_opt = 0 ! meant to redistribute the points in z
@@ -55,6 +64,10 @@ implicit none
REAL(dp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition
INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R
+ ! pb_phase, for parallel boundary phase, contains the factor that occurs when taking into account
+ ! that q0 is defined in the middle of the fluxtube whereas the radial position spans in [0,Lx)
+ ! This shift introduces a (-1)^(Nexc*iky) phase change that is included in GENE
+ COMPLEX(dp), PUBLIC, DIMENSION(:), ALLOCATABLE :: pb_phase_L, pb_phase_R
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
@@ -66,9 +79,19 @@ CONTAINS
! Read the input parameters
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
- kappa, s_kappa,delta, s_delta, zeta, s_zeta ! For miller
+ kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For miller
+ parallel_bc, shift_y
READ(lu_in,geometry)
IF(shear .NE. 0._dp) SHEARED = .true.
+ SELECT CASE(parallel_bc)
+ CASE ('dirichlet')
+ CASE ('periodic')
+ CASE ('shearless')
+ CASE ('disconnected')
+ CASE DEFAULT
+ stop 'Parallel BC not recognized'
+ END SELECT
+ IF(my_id .EQ. 0) print*, 'Parallel BC : ', parallel_bc
END SUBROUTINE geometry_readinputs
@@ -324,12 +347,13 @@ CONTAINS
SUBROUTINE set_ikx_zBC_map
IMPLICIT NONE
REAL :: shift, kx_shift
- ! For periodic CHI BC or 0 dirichlet
- LOGICAL :: PERIODIC_CHI_BC = .FALSE.
- ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
- ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
+ INTEGER :: ikx_shift
- !! No shear case (simple id mapping)
+ ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
+ ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
+ ALLOCATE(pb_phase_L(ikys:ikye))
+ ALLOCATE(pb_phase_R(ikys:ikye))
+ !! No shear case (simple id mapping) or not at the end of the z domain
!3 | 1 2 3 4 5 6 | ky = 3 dky
!2 ky | 1 2 3 4 5 6 | ky = 2 dky
!1 A | 1 2 3 4 5 6 | ky = 1 dky
@@ -337,65 +361,109 @@ CONTAINS
!(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=free)
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
- ikx_zBC_L(iky,ikx) = ikx
+ ikx_zBC_L(iky,ikx) = ikx ! connect to itself per default
ikx_zBC_R(iky,ikx) = ikx
ENDDO
+ pb_phase_L(iky) = 1._dp ! no phase change per default
+ pb_phase_R(iky) = 1._dp
ENDDO
- ! Modify connection map only at border of z
- IF(SHEARED) THEN
- ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (even NZ)
- !3 | 4 x x x 2 3 | ky = 3 dky
- !2 ky | 3 4 x x 1 2 | ky = 2 dky
- !1 A | 2 3 4 x 6 1 | ky = 1 dky
- !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
- !kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
-
- ! connection map BC of the RIGHT boundary (z=pi*Npol-dz) (ODD NZ)
- !3 | x x x 2 3 | ky = 3 dky
- !2 ky | 3 x x 1 2 | ky = 2 dky
- !1 A | 2 3 x 5 1 | ky = 1 dky
- !0 | -> kx | 1____2____3____4____5 | ky = 0 dky
- !kx = 0 0.1 0.2 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
- IF(contains_zmax) THEN ! Check if the process is at the end of the FT
+ ! Parallel boundary are not trivial for sheared case and if
+ ! the user does not ask explicitly for shearless bc
+ IF(SHEARED .AND. (parallel_bc .NE. 'shearless')) THEN
+ !!!!!!!!!! LEFT PARALLEL BOUNDARY
+ ! 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 = ikys,ikye
+ ! Formula for the shift due to shear after Npol turns
shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
- kx_shift = kxarray(ikx) + shift
- ! We use EPSILON() to treat perfect equality case
- IF( ((kx_shift-EPSILON(kx_shift)) .GT. kx_max) .AND. (.NOT. PERIODIC_CHI_BC) )THEN ! outside of the frequ domain
- ikx_zBC_R(iky,ikx) = -99
- ELSE
- ikx_zBC_R(iky,ikx) = ikx+(iky-1)*Nexc
- IF( ikx_zBC_R(iky,ikx) .GT. Nkx ) &
- ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - 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') !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)
ENDDO
ENDIF
- ! connection map BC of the LEFT boundary (z=-pi*Npol)
- !3 | x 5 6 1 x x | ky = 3 dky
- !2 ky | 5 6 1 2 x x | ky = 2 dky
- !1 A | 6 1 2 3 x 5 | ky = 1 dky
- !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=2pi*shear*npol*dky)
- IF(contains_zmin) THEN ! Check if the process is at the start of the FT
+ ! Option for disconnecting every modes, viz. connecting all boundary to 0
+ IF(parallel_bc .EQ. 'disconnected') ikx_zBC_L = -99
+ !!!!!!!!!! RIGHT PARALLEL BOUNDARY
+ IF(contains_zmax) THEN ! Check if the process is at the end of the flux-tube
DO iky = ikys,ikye
+ ! Formula for the shift due to shear after Npol
shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
- kx_shift = kxarray(ikx) - shift
- ! We use EPSILON() to treat perfect equality case
- IF( ((kx_shift+EPSILON(kx_shift)) .LT. kx_min) .AND. (.NOT. PERIODIC_CHI_BC) ) THEN ! outside of the frequ domain
- ikx_zBC_L(iky,ikx) = -99
- ELSE
- ikx_zBC_L(iky,ikx) = ikx-(iky-1)*Nexc
- IF( ikx_zBC_L(iky,ikx) .LT. 1 ) &
- ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + 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') !reroute it by cycling through modes
+ 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)
ENDDO
ENDIF
- ELSE
-ENDIF
+ ! Option for disconnecting every modes, viz. connecting all boundary to 0
+ IF(parallel_bc .EQ. 'disconnected') ikx_zBC_R = -99
+ ENDIF
+ ! write(*,*) kxarray
+ ! write(*,*) kyarray
+ ! write(*,*) 'ikx_zBC_L :-----------'
+ ! DO iky = ikys,ikye
+ ! print*, ikx_zBC_L(iky,:)
+ ! enddo
+ ! write(*,*) 'ikx_zBC_R :-----------'
+ ! DO iky = ikys,ikye
+ ! print*, ikx_zBC_R(iky,:)
+ ! enddo
+ ! stop
+ !!!!!!! Example of maps ('x' means connected to 0 value, in the table it is -99)
+ ! dirichlet connection map BC of the RIGHT boundary (z=pi*Npol-dz)
+ !3 | 4 x x x 2 3 | ky = 3 dky
+ !2 ky | 3 4 x x 1 2 | ky = 2 dky
+ !1 A | 2 3 4 x 6 1 | ky = 1 dky
+ !0 | -> kx | 1____2____3____4____5____6 | ky = 0 dky
+ !kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
+
+ ! periodic connection map BC of the LEFT boundary (z=-pi*Npol)
+ !3 | 4 5 6 1 2 3 | ky = 3 dky
+ !2 ky | 5 6 1 2 3 4 | ky = 2 dky
+ !1 A | 6 1 2 3 4 5 | ky = 1 dky
+ !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=2pi*shear*npol*dky)
+
+ ! shearless connection map BC of the LEFT/RIGHT boundary (z=+/-pi*Npol)
+ !3 | 1 2 3 4 5 6 | ky = 3 dky
+ !2 ky | 1 2 3 4 5 6 | ky = 2 dky
+ !1 A | 1 2 3 4 5 6 | ky = 1 dky
+ !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=2pi*shear*npol*dky)
+
+ ! disconnected connection map BC of the LEFT/RIGHT boundary (z=+/-pi*Npol)
+ !3 | x x x x x x | ky = 3 dky
+ !2 ky | x x x x x x | ky = 2 dky
+ !1 A | x x x x x x | ky = 1 dky
+ !0 | -> kx | x____x____x____x____x____x | ky = 0 dky
+ !(e.g.) kx = 0 0.1 0.2 0.3 -0.2 -0.1 (dkx=2pi*shear*npol*dky)
END SUBROUTINE set_ikx_zBC_map
!
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 0f3136bf665f2be95b9f90c265859ad4cf40d990..bd23ad9217f795eefb7489e3b37387627c4fc616 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -3,7 +3,7 @@ USE basic
USE grid
USE time_integration
USE model, ONLY : KIN_E, beta
-USE geometry, ONLY : SHEARED, ikx_zBC_L, ikx_zBC_R
+USE geometry, ONLY : SHEARED, ikx_zBC_L, ikx_zBC_R, pb_phase_L, pb_phase_R
IMPLICIT NONE
INTEGER :: status(MPI_STATUS_SIZE), source, dest, count, ipg
@@ -39,7 +39,7 @@ SUBROUTINE update_ghosts_EM
IF(beta .GT. 0._dp) &
CALL update_ghosts_z_psi
ENDIF
-
+
CALL cpu_time(t1_ghost)
tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
END SUBROUTINE update_ghosts_EM
@@ -159,9 +159,9 @@ SUBROUTINE update_ghosts_z_e
ikxBC_L = ikx_zBC_L(iky,ikx);
IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! first-1 gets last
- moments_e(:,:,iky,ikx,izs-1,updatetlevel) = buff_pjxy_zBC_e(:,:,iky,ikxBC_L,-1)
+ moments_e(:,:,iky,ikx,izs-1,updatetlevel) = pb_phase_L(iky)*buff_pjxy_zBC_e(:,:,iky,ikxBC_L,-1)
! first-2 gets last-1
- moments_e(:,:,iky,ikx,izs-2,updatetlevel) = buff_pjxy_zBC_e(:,:,iky,ikxBC_L,-2)
+ moments_e(:,:,iky,ikx,izs-2,updatetlevel) = pb_phase_L(iky)*buff_pjxy_zBC_e(:,:,iky,ikxBC_L,-2)
ELSE
moments_e(:,:,iky,ikx,izs-1,updatetlevel) = 0._dp
moments_e(:,:,iky,ikx,izs-2,updatetlevel) = 0._dp
@@ -169,9 +169,9 @@ SUBROUTINE update_ghosts_z_e
ikxBC_R = ikx_zBC_R(iky,ikx);
IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
- moments_e(:,:,iky,ikx,ize+1,updatetlevel) = buff_pjxy_zBC_e(:,:,iky,ikxBC_R,+1)
+ moments_e(:,:,iky,ikx,ize+1,updatetlevel) = pb_phase_R(iky)*buff_pjxy_zBC_e(:,:,iky,ikxBC_R,+1)
! last+2 gets first+1
- moments_e(:,:,iky,ikx,ize+2,updatetlevel) = buff_pjxy_zBC_e(:,:,iky,ikxBC_R,+2)
+ moments_e(:,:,iky,ikx,ize+2,updatetlevel) = pb_phase_R(iky)*buff_pjxy_zBC_e(:,:,iky,ikxBC_R,+2)
ELSE
moments_e(:,:,iky,ikx,ize+1,updatetlevel) = 0._dp
moments_e(:,:,iky,ikx,ize+2,updatetlevel) = 0._dp
@@ -217,9 +217,9 @@ SUBROUTINE update_ghosts_z_i
ikxBC_L = ikx_zBC_L(iky,ikx);
IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! first-1 gets last
- moments_i(:,:,iky,ikx,izs-1,updatetlevel) = buff_pjxy_zBC_i(:,:,iky,ikxBC_L,-1)
+ moments_i(:,:,iky,ikx,izs-1,updatetlevel) = pb_phase_L(iky)*buff_pjxy_zBC_i(:,:,iky,ikxBC_L,-1)
! first-2 gets last-1
- moments_i(:,:,iky,ikx,izs-2,updatetlevel) = buff_pjxy_zBC_i(:,:,iky,ikxBC_L,-2)
+ moments_i(:,:,iky,ikx,izs-2,updatetlevel) = pb_phase_L(iky)*buff_pjxy_zBC_i(:,:,iky,ikxBC_L,-2)
ELSE
moments_i(:,:,iky,ikx,izs-1,updatetlevel) = 0._dp
moments_i(:,:,iky,ikx,izs-2,updatetlevel) = 0._dp
@@ -227,9 +227,9 @@ SUBROUTINE update_ghosts_z_i
ikxBC_R = ikx_zBC_R(iky,ikx);
IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
- moments_i(:,:,iky,ikx,ize+1,updatetlevel) = buff_pjxy_zBC_i(:,:,iky,ikxBC_R,+1)
+ moments_i(:,:,iky,ikx,ize+1,updatetlevel) = pb_phase_R(iky)*buff_pjxy_zBC_i(:,:,iky,ikxBC_R,+1)
! last+2 gets first+1
- moments_i(:,:,iky,ikx,ize+2,updatetlevel) = buff_pjxy_zBC_i(:,:,iky,ikxBC_R,+2)
+ moments_i(:,:,iky,ikx,ize+2,updatetlevel) = pb_phase_R(iky)*buff_pjxy_zBC_i(:,:,iky,ikxBC_R,+2)
ELSE
moments_i(:,:,iky,ikx,ize+1,updatetlevel) = 0._dp
moments_i(:,:,iky,ikx,ize+2,updatetlevel) = 0._dp
@@ -277,9 +277,9 @@ SUBROUTINE update_ghosts_z_phi
ikxBC_L = ikx_zBC_L(iky,ikx);
IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! first-1 gets last
- phi(iky,ikx,izs-1) = buff_xy_zBC(iky,ikxBC_L,-1)
+ phi(iky,ikx,izs-1) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-1)
! first-2 gets last-1
- phi(iky,ikx,izs-2) = buff_xy_zBC(iky,ikxBC_L,-2)
+ phi(iky,ikx,izs-2) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-2)
ELSE
phi(iky,ikx,izs-1) = 0._dp
phi(iky,ikx,izs-2) = 0._dp
@@ -287,9 +287,9 @@ SUBROUTINE update_ghosts_z_phi
ikxBC_R = ikx_zBC_R(iky,ikx);
IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
- phi(iky,ikx,ize+1) = buff_xy_zBC(iky,ikxBC_R,+1)
+ phi(iky,ikx,ize+1) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,+1)
! last+2 gets first+1
- phi(iky,ikx,ize+2) = buff_xy_zBC(iky,ikxBC_R,+2)
+ phi(iky,ikx,ize+2) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,+2)
ELSE
phi(iky,ikx,ize+1) = 0._dp
phi(iky,ikx,ize+2) = 0._dp
@@ -339,9 +339,9 @@ SUBROUTINE update_ghosts_z_psi
ikxBC_L = ikx_zBC_L(iky,ikx);
IF (ikxBC_L .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! first-1 gets last
- psi(iky,ikx,izs-1) = buff_xy_zBC(iky,ikxBC_L,-1)
+ psi(iky,ikx,izs-1) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-1)
! first-2 gets last-1
- psi(iky,ikx,izs-2) = buff_xy_zBC(iky,ikxBC_L,-2)
+ psi(iky,ikx,izs-2) = pb_phase_L(iky)*buff_xy_zBC(iky,ikxBC_L,-2)
ELSE
psi(iky,ikx,izs-1) = 0._dp
psi(iky,ikx,izs-2) = 0._dp
@@ -349,9 +349,9 @@ SUBROUTINE update_ghosts_z_psi
ikxBC_R = ikx_zBC_R(iky,ikx);
IF (ikxBC_R .NE. -99) THEN ! Exchanging the modes that have a periodic pair (a)
! last+1 gets first
- psi(iky,ikx,ize+1) = buff_xy_zBC(iky,ikxBC_R,+1)
+ psi(iky,ikx,ize+1) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,+1)
! last+2 gets first+1
- psi(iky,ikx,ize+2) = buff_xy_zBC(iky,ikxBC_R,+2)
+ psi(iky,ikx,ize+2) = pb_phase_R(iky)*buff_xy_zBC(iky,ikxBC_R,+2)
ELSE
psi(iky,ikx,ize+1) = 0._dp
psi(iky,ikx,ize+2) = 0._dp