From de0ac4c5f642f8af21f66d8bec8132e443a03a3b Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Thu, 23 Jun 2022 16:11:33 +0200
Subject: [PATCH] first shear implementation
---
src/geometry_mod.F90 | 64 +++++++++++++++++++++--
src/ghosts_mod.F90 | 122 ++++++++++++++++++++++++++++++++++++-------
2 files changed, 162 insertions(+), 24 deletions(-)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index bd67f738..ad36cc5c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -18,7 +18,7 @@ implicit none
REAL(dp), PUBLIC, PROTECTED :: q0 = 1.4_dp ! safety factor
REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
REAL(dp), PUBLIC, PROTECTED :: eps = 0.18_dp ! inverse aspect ratio
-
+ LOGICAL, PUBLIC, PROTECTED :: SHEARED = .false. ! flag for shear magn. geom or not
! Geometrical operators
! Curvature
REAL(dp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
@@ -36,10 +36,11 @@ implicit none
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
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_map
! Functions
- PUBLIC :: geometry_readinputs, geometry_outputinputs, eval_magnetic_geometry
-
+ PUBLIC :: geometry_readinputs, geometry_outputinputs,&
+ eval_magnetic_geometry, set_ikx_zBC_map
CONTAINS
@@ -48,6 +49,7 @@ CONTAINS
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps
READ(lu_in,geometry)
+ IF(shear .NE. 0._dp) SHEARED = .true.
END SUBROUTINE geometry_readinputs
@@ -93,6 +95,8 @@ CONTAINS
ENDDO
ENDDO
ENDDO
+ ! set the mapping for parallel boundary conditions
+ CALL set_ikx_zBC_map
two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
!
@@ -259,10 +263,62 @@ CONTAINS
ENDDO parity
END SUBROUTINE eval_1D_geometry
+
!
!--------------------------------------------------------------------------------
!
+ SUBROUTINE set_ikx_zBC_map
+ IMPLICIT NONE
+ INTEGER, DIMENSION(:), ALLOCATABLE :: ikx_array
+ INTEGER :: shift
+ ALLOCATE(ikx_array(ikxs:ikxe))
+ DO ikx = ikxs,ikxe
+ ikx_array(ikx) = MODULO(ikx - Nx/2,Nx) + 1
+ ENDDO
+ IF(SHEARED) THEN
+ !! We allocate a mapping to tell where the current mode will point for the
+ ! parallel periodic sheared BC for the kx index:
+ ! map for 0 Dirichlet BC
+ !3 | 1 2 3 4 5 -1 -1 -1|
+ !2 ky | 8 1 2 3 4 5 -1 -1|
+ !1 A | 7 8 1 2 3 4 5 -1|
+ !0 | -> kx | 6____7____8____1____2____3____4____5|
+ ! map for periodic BC
+ !3 | 1 2 3 4 5 -1 -1 -1|
+ !2 ky | 8 1 2 3 4 5 -1 -1|
+ !1 A | 7 8 1 2 3 4 5 -1|
+ !0 | -> kx | 6____7____8____1____2____3____4____5|
+ ALLOCATE(ikx_zBC_map(ikys:ikye,ikxs:ikxe))
+ ikx_zBC_map(ikys:ikye,ikxs:ikxe) = -1
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe - iky + 1
+ shift = ikx_array(MODULO(ikx+iky-1,Nx+1))
+ ikx_zBC_map(iky,ikx) = shift
+ ENDDO
+ ENDDO
+ ELSE
+ !! No shear case (simple id mapping)
+ !3 | 6 7 8 1 2 3 4 5|
+ !2 ky | 6 7 8 1 2 3 4 5|
+ !1 A | 6 7 8 1 2 3 4 5|
+ !0 | -> kx | 6____7____8____1____2____3____4____5|
+ DO iky = ikys,ikye
+ ikx_zBC_map(iky,:) = ikx_array(:)
+ ENDDO
+ENDIF
+! IF (my_id .EQ. 0) THEN
+! write(*,*) 'ikx map for parallel BC'
+! DO, iky = ikys,ikye
+! write(*,*) ( ikx_zBC_map(iky,ikx), ikx=ikxs,ikxe )
+! enddo
+! ENDIF
+END SUBROUTINE set_ikx_zBC_map
+
+!
+!--------------------------------------------------------------------------------
+!
+
SUBROUTINE geometry_allocate_mem
! Curvature and geometry
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 9dff68f4..14da9153 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -4,6 +4,7 @@ USE fields, ONLY : moments_e, moments_i, phi
USE grid
USE time_integration
USE model, ONLY : KIN_E
+USE geometry, ONLY : SHEARED, ikx_zBC_map
IMPLICIT NONE
@@ -122,7 +123,7 @@ END SUBROUTINE update_ghosts_z_moments
SUBROUTINE update_ghosts_z_e
IMPLICIT NONE
-
+ INTEGER :: ikxBC
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF (num_procs_z .GT. 1) THEN
@@ -145,20 +146,44 @@ SUBROUTINE update_ghosts_z_e
moments_e(:,:,:,:,ize+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 23, & ! Recieve from Up the last+2
comm0, status, ierr)
ELSE ! still need to perform periodic boundary conditions
- ! first-1 gets last
- moments_e(:,:,:,:,izs-1,updatetlevel) = moments_e(:,:,:,:,ize ,updatetlevel)
- ! first-2 gets last-1
- moments_e(:,:,:,:,izs-2,updatetlevel) = moments_e(:,:,:,:,ize-1,updatetlevel)
- ! last+1 gets first
- moments_e(:,:,:,:,ize+1,updatetlevel) = moments_e(:,:,:,:,izs ,updatetlevel)
- ! last+2 gets first+1
- moments_e(:,:,:,:,ize+2,updatetlevel) = moments_e(:,:,:,:,izs+1,updatetlevel)
+ IF(SHEARED) THEN
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ ikxBC = ikx_zBC_map(ikx,iky);
+ IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a)
+ ! first-1 gets last
+ moments_e(iky,ikx,:,:,izs-1,updatetlevel) = moments_e(iky,ikxBC,:,:,ize ,updatetlevel)
+ ! first-2 gets last-1
+ moments_e(iky,ikx,:,:,izs-2,updatetlevel) = moments_e(iky,ikxBC,:,:,ize-1,updatetlevel)
+ ! last+1 gets first
+ moments_e(iky,ikx,:,:,ize+1,updatetlevel) = moments_e(iky,ikxBC,:,:,izs ,updatetlevel)
+ ! last+2 gets first+1
+ moments_e(iky,ikx,:,:,ize+2,updatetlevel) = moments_e(iky,ikxBC,:,:,izs+1,updatetlevel)
+ ELSE
+ moments_e(iky,ikx,:,:,izs-1,updatetlevel) = 0._dp
+ moments_e(iky,ikx,:,:,izs-2,updatetlevel) = 0._dp
+ moments_e(iky,ikx,:,:,ize+1,updatetlevel) = 0._dp
+ moments_e(iky,ikx,:,:,ize+2,updatetlevel) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ ELSE ! No shear so simple periodic BC
+ ! first-1 gets last
+ moments_e(:,:,:,:,izs-1,updatetlevel) = moments_e(:,:,:,:,ize ,updatetlevel)
+ ! first-2 gets last-1
+ moments_e(:,:,:,:,izs-2,updatetlevel) = moments_e(:,:,:,:,ize-1,updatetlevel)
+ ! last+1 gets first
+ moments_e(:,:,:,:,ize+1,updatetlevel) = moments_e(:,:,:,:,izs ,updatetlevel)
+ ! last+2 gets first+1
+ moments_e(:,:,:,:,ize+2,updatetlevel) = moments_e(:,:,:,:,izs+1,updatetlevel)
+ ENDIF
ENDIF
END SUBROUTINE update_ghosts_z_e
SUBROUTINE update_ghosts_z_i
IMPLICIT NONE
+ INTEGER :: ikxBC
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF (num_procs_z .GT. 1) THEN
@@ -180,20 +205,45 @@ SUBROUTINE update_ghosts_z_i
CALL mpi_sendrecv(moments_i(:,:,:,:,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 33, & ! Send to Down the first+1
moments_i(:,:,:,:,ize+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 33, & ! Recieve from Down the last+2
comm0, status, ierr)
- ELSE ! still need to perform periodic boundary conditions
- moments_i(:,:,:,:,izs-1,updatetlevel) = moments_i(:,:,:,:,ize ,updatetlevel)
-
- moments_i(:,:,:,:,izs-2,updatetlevel) = moments_i(:,:,:,:,ize-1,updatetlevel)
-
- moments_i(:,:,:,:,ize+1,updatetlevel) = moments_i(:,:,:,:,izs ,updatetlevel)
-
- moments_i(:,:,:,:,ize+2,updatetlevel) = moments_i(:,:,:,:,izs+1,updatetlevel)
- ENDIF
+ ELSE ! still need to perform periodic boundary conditions
+ IF(SHEARED) THEN
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ ikxBC = ikx_zBC_map(iky,ikx);
+ IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a)
+ ! first-1 gets last
+ moments_i(:,:,iky,ikx,izs-1,updatetlevel) = moments_i(:,:,iky,ikxBC,ize ,updatetlevel)
+ ! first-2 gets last-1
+ moments_i(:,:,iky,ikx,izs-2,updatetlevel) = moments_i(:,:,iky,ikxBC,ize-1,updatetlevel)
+ ! last+1 gets first
+ moments_i(:,:,iky,ikx,ize+1,updatetlevel) = moments_i(:,:,iky,ikxBC,izs ,updatetlevel)
+ ! last+2 gets first+1
+ moments_i(:,:,iky,ikx,ize+2,updatetlevel) = moments_i(:,:,iky,ikxBC,izs+1,updatetlevel)
+ ELSE
+ moments_i(:,:,iky,ikx,izs-1,updatetlevel) = 0._dp
+ moments_i(:,:,iky,ikx,izs-2,updatetlevel) = 0._dp
+ moments_i(:,:,iky,ikx,ize+1,updatetlevel) = 0._dp
+ moments_i(:,:,iky,ikx,ize+2,updatetlevel) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ ELSE ! No shear so simple periodic BC
+ ! first-1 gets last
+ moments_i(:,:,:,:,izs-1,updatetlevel) = moments_i(:,:,:,:,ize ,updatetlevel)
+ ! first-2 gets last-1
+ moments_i(:,:,:,:,izs-2,updatetlevel) = moments_i(:,:,:,:,ize-1,updatetlevel)
+ ! last+1 gets first
+ moments_i(:,:,:,:,ize+1,updatetlevel) = moments_i(:,:,:,:,izs ,updatetlevel)
+ ! last+2 gets first+1
+ moments_i(:,:,:,:,ize+2,updatetlevel) = moments_i(:,:,:,:,izs+1,updatetlevel)
+ ENDIF
+ ENDIF
END SUBROUTINE update_ghosts_z_i
SUBROUTINE update_ghosts_z_phi
- IMPLICIT NONE
+ IMPLICIT NONE
+ INTEGER :: ikxBC
CALL cpu_time(t1_ghost)
IF(Nz .GT. 1) THEN
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
@@ -220,7 +270,39 @@ SUBROUTINE update_ghosts_z_phi
phi(:,:,izs-1) = phi(:,:,ize)
phi(:,:,izs-2) = phi(:,:,ize-1)
phi(:,:,ize+1) = phi(:,:,izs)
- phi(:,:,ize+2) = phi(:,:,izs+1)
+ phi(:,:,ize+2) = phi(:,:,izs+1)
+
+ IF(SHEARED) THEN
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ ikxBC = ikx_zBC_map(iky,ikx);
+ IF (ikxBC .NE. -1) THEN ! Exchanging the modes that have a periodic pair (a)
+ ! first-1 gets last
+ phi(iky,ikx,izs-1) = phi(iky,ikxBC,ize )
+ ! first-2 gets last-1
+ phi(iky,ikx,izs-2) = phi(iky,ikxBC,ize-1)
+ ! last+1 gets first
+ phi(iky,ikx,ize+1) = phi(iky,ikxBC,izs )
+ ! last+2 gets first+1
+ phi(iky,ikx,ize+2) = phi(iky,ikxBC,izs+1)
+ ELSE
+ phi(iky,ikx,izs-1) = 0._dp
+ phi(iky,ikx,izs-2) = 0._dp
+ phi(iky,ikx,ize+1) = 0._dp
+ phi(iky,ikx,ize+2) = 0._dp
+ ENDIF
+ ENDDO
+ ENDDO
+ ELSE ! No shear so simple periodic BC
+ ! first-1 gets last
+ phi(:,:,izs-1) = phi(:,:,ize )
+ ! first-2 gets last-1
+ phi(:,:,izs-2) = phi(:,:,ize-1)
+ ! last+1 gets first
+ phi(:,:,ize+1) = phi(:,:,izs )
+ ! last+2 gets first+1
+ phi(:,:,ize+2) = phi(:,:,izs+1)
+ ENDIF
ENDIF
ENDIF
CALL cpu_time(t1_ghost)
--
GitLab