From 37641e0d94b6dfc1918e2d3b2a4353c02d3a6fa8 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 21 Sep 2022 11:18:05 +0200
Subject: [PATCH] Clean and Added Nexc parameter, not working yet for Nexc > 1
---
src/geometry_mod.F90 | 60 +++++++++++++++++++++++---------------------
1 file changed, 32 insertions(+), 28 deletions(-)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index ff0e309e..74ec32c5 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -151,25 +151,18 @@ CONTAINS
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = eps * SIN(z) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
-
- ! ! Curvature operator
- ! DO iky = ikys, ikye
- ! ky = kyarray(iky)
- ! DO ikx= ikxs, ikxe
- ! kx = kxarray(ikx)
- ! Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - ( COS(z) + (shear*z-alpha_MHD*SIN(z))*SIN(z) )*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- ! ENDDO
- ! ENDDO
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = eps*SIN(z)/hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
+
! Curvature operator
Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Kx
- Gy = -COS(Z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky
+ Gy = -COS(z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)* hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ ! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
@@ -446,9 +439,15 @@ CONTAINS
SUBROUTINE set_ikx_zBC_map
IMPLICIT NONE
- REAL :: shift, kx_shift, kxmax_, kxmin_
+ 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))
+
+ ! No periodic connection for extension of the domain
+ IF(Nexc .GT. 1) PERIODIC_CHI_BC = .false.
+
!! No shear case (simple id mapping)
!3 | 1 2 3 4 5 6 | ky = 3 dky
!2 ky | 1 2 3 4 5 6 | ky = 2 dky
@@ -463,25 +462,31 @@ CONTAINS
ENDDO
! Modify connection map only at border of z
IF(SHEARED) THEN
- ! connection map BC of the RIGHT boundary (z=pi*Npol-dz)
+ ! 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)
- kxmax_ = kx_max
+
+ ! 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
DO iky = ikys,ikye
- shift = 2._dp*PI*shear*kyarray(iky)*Npol*Nexc
+ shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) + shift
- IF(kx_shift .GT. kxmax_) THEN ! outside of the frequ domain
+ ! 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
- IF(SINGLE_KY) ikx_zBC_R(iky,ikx) = (ikx+(iky+1))-1
+ 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
+ ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - Nkx
ENDIF
ENDDO
ENDDO
@@ -492,19 +497,18 @@ CONTAINS
!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)
- kxmin_ = -kx_max+deltakx
IF(contains_zmin) THEN ! Check if the process is at the start of the FT
DO iky = ikys,ikye
- shift = 2._dp*PI*shear*kyarray(iky)*Npol*Nexc
+ shift = 2._dp*PI*shear*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) - shift
- IF( kx_shift .LT. kxmin_ ) THEN ! outside of the frequ domain
+ ! 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
- IF(SINGLE_KY) ikx_zBC_L(iky,ikx) = (ikx-(iky+1))+1
- IF( ikx_zBC_L(iky,ikx) .LE. 0 ) &
- ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + Nkx
+ 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
ENDIF
ENDDO
ENDDO
--
GitLab