module geometry
! computes geometrical quantities
! Adapted from B.J.Frei MOLIX code (2021)
use prec_const
use model
use grid
use array
use fields
use basic
use calculus, ONLY: simpson_rule_z
implicit none
PRIVATE
! Geometry parameters
CHARACTER(len=16), &
PUBLIC, PROTECTED :: geom
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
! Jacobian
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
COMPLEX(dp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
! Metric
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
! Relative magnetic field strength
REAL(dp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, hatB_NL
! Relative strength of major radius
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_L, ikx_zBC_R
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
eval_magnetic_geometry, set_ikx_zBC_map
CONTAINS
SUBROUTINE geometry_readinputs
! Read the input parameters
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps
READ(lu_in,geometry)
IF(shear .NE. 0._dp) SHEARED = .true.
END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
COMPLEX(dp), DIMENSION(izs:ize) :: integrant
INTEGER :: fid
! Allocate arrays
CALL geometry_allocate_mem
!
IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
call eval_1D_geometry
ELSE
SELECT CASE(geom)
CASE('circular')
IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
call eval_circular_geometry
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
CASE('Z-pinch')
IF( my_id .eq. 0 ) WRITE(*,*) 'Z-pinch geometry'
call eval_zpinch_geometry
SHEARED = .FALSE.
CASE DEFAULT
ERROR STOP 'Error stop: geometry not recognized!!'
END SELECT
ENDIF
!
! Evaluate perpendicular wavenumber
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
DO eo = 0,1
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx = ikxs, ikxe
kx = kxarray(ikx)
DO iz = izgs,izge
kparray(iky, ikx, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
ENDDO
ENDDO
ENDDO
ENDDO
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
integrant = Jacobian(izs:ize,0) ! Convert into complex array
CALL simpson_rule_z(integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z, kx, ky, alpha_MHD
alpha_MHD = 0._dp
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = shear*z - alpha_MHD*SIN(z)
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) !avoid 1/0 in Zpinch config
gzz(iz,eo) = 0._dp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
hatR(iz,eo) = 1._dp + eps*COS(z)
hatZ(iz,eo) = 1._dp + eps*SIN(z)
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
hatB (iz,eo) = 1._dp / hatR(iz,eo)
hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
! 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
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_salphaB_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_circular_geometry
! evaluate circular geometry model
! Ref: Lapilonne et al., PoP, 2009
implicit none
REAL(dp) :: X, kx, ky, Gamma1, Gamma2, Gamma3
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
X = zarray(iz,eo) - eps*SIN(zarray(iz,eo)) ! chi = theta - eps sin(theta)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = shear*X - eps*SIN(X)
gxz(iz,eo) = - SIN(X)
gyy(iz,eo) = 1._dp + (shear*X)**2 - 2._dp*eps*COS(X) - 2._dp*shear*X*eps*SIN(X)
gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) - 2._dp*COS(X) - shear*X*SIN(X)
gzz(iz,eo) = 1._dp/eps**2 - 2._dp*COS(X)/eps
dxdR(iz,eo)= COS(X)
dxdZ(iz,eo)= SIN(X)
! Relative strengh of radius
hatR(iz,eo) = 1._dp + eps*COS(X)
hatZ(iz,eo) = 1._dp + eps*SIN(X)
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = X
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)**2
! Relative strengh of modulus of B
hatB (iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2)
hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -(COS(X) + eps*SIN(X)**2)/hatB(iz,eo)**2 ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = eps * SIN(X) * (1._dp - eps*COS(X)) / hatB(iz,eo)**2 ! Gene put a factor hatB or 1/hatR in this
Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo)
Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)
Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(iky, ikx, iz,eo) = Gamma1*((kx + shear*X*ky)*gradzB(iz,eo)/eps - gradxB(iz,eo)*ky) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_circular_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE eval_zpinch_geometry
! evaluate s-alpha geometry model
implicit none
REAL(dp) :: z, kx, ky, alpha_MHD
alpha_MHD = 0._dp
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = 0._dp
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp
gyz(iz,eo) = 0._dp
gzz(iz,eo) = 1._dp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
hatR(iz,eo) = 1._dp
hatZ(iz,eo) = 1._dp
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
hatB (iz,eo) = 1._dp
hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = 0._dp ! 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) = -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
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_zpinch_geometry
!
!--------------------------------------------------------------------------------
!
subroutine eval_1D_geometry
! evaluate 1D perp geometry model
implicit none
REAL(dp) :: z, kx, ky
parity: DO eo = 0,1
zloop: DO iz = izs,ize
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
gxy(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp
! Relative strengh of radius
hatR(iz,eo) = 1._dp
! Jacobian
Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
hatB(iz,eo) = 1._dp
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(ikx, iky, iz,eo) = -kx ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp
ENDDO zloop
ENDDO parity
END SUBROUTINE eval_1D_geometry
!
!--------------------------------------------------------------------------------
!
SUBROUTINE set_ikx_zBC_map
IMPLICIT NONE
REAL :: shift, kx_shift, kxmax_, kxmin_
ALLOCATE(ikx_zBC_R(ikys:ikye,ikxs:ikxe))
ALLOCATE(ikx_zBC_L(ikys:ikye,ikxs:ikxe))
!! 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
!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=free)
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
ikx_zBC_L(iky,ikx) = ikx
ikx_zBC_R(iky,ikx) = ikx
ENDDO
ENDDO
! Modify connection map only at border of z
IF(SHEARED) THEN
! 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)
kxmax_ = kx_max
IF(contains_zmax) THEN ! Check if the process is at the end of the FT
DO iky = ikys,ikye
shift = 2._dp*shear*PI*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) + shift
IF(kx_shift .GT. kxmax_) 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
IF( ikx_zBC_R(iky,ikx) .GT. Nkx ) &
ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - Nkx
ENDIF
ENDDO
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)
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*shear*PI*kyarray(iky)*Npol
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) - shift
IF( kx_shift .LT. kxmin_ ) 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
ENDIF
ENDDO
ENDDO
ENDIF
ELSE
ENDIF
END SUBROUTINE set_ikx_zBC_map
!
!--------------------------------------------------------------------------------
!
SUBROUTINE geometry_allocate_mem
! Curvature and geometry
CALL allocate_array( Ckxky, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
CALL allocate_array( Jacobian,izgs,izge, 0,1)
CALL allocate_array( gxx,izgs,izge, 0,1)
CALL allocate_array( gxy,izgs,izge, 0,1)
CALL allocate_array( gxz,izgs,izge, 0,1)
CALL allocate_array( gyy,izgs,izge, 0,1)
CALL allocate_array( gyz,izgs,izge, 0,1)
CALL allocate_array( gzz,izgs,izge, 0,1)
CALL allocate_array( gradxB,izgs,izge, 0,1)
CALL allocate_array( gradyB,izgs,izge, 0,1)
CALL allocate_array( gradzB,izgs,izge, 0,1)
CALL allocate_array( hatB,izgs,izge, 0,1)
CALL allocate_array( hatB_NL,izgs,izge, 0,1)
CALL allocate_array( hatR,izgs,izge, 0,1)
CALL allocate_array( hatZ,izgs,izge, 0,1)
CALL allocate_array( Rc,izgs,izge, 0,1)
CALL allocate_array( phic,izgs,izge, 0,1)
CALL allocate_array( Zc,izgs,izge, 0,1)
CALL allocate_array( dxdR,izgs,izge, 0,1)
CALL allocate_array( dxdZ,izgs,izge, 0,1)
call allocate_array(gradz_coeff,izgs,izge, 0,1)
CALL allocate_array( kparray, ikys,ikye, ikxs,ikxe,izgs,izge,0,1)
END SUBROUTINE geometry_allocate_mem
SUBROUTINE geometry_outputinputs(fidres, str)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach
USE prec_const
IMPLICIT NONE
INTEGER, INTENT(in) :: fidres
CHARACTER(len=256), INTENT(in) :: str
CALL attach(fidres, TRIM(str),"geometry", geom)
CALL attach(fidres, TRIM(str), "q0", q0)
CALL attach(fidres, TRIM(str), "shear", shear)
CALL attach(fidres, TRIM(str), "eps", eps)
END SUBROUTINE geometry_outputinputs
end module geometry