Homec4science
Diffusion GYACOMO (OBSOLETE: LATEST VERSION ON gitlab.epfl.ch) (half_ky_plane) Browse

gyacomo/src/geometry_mod.F90
8dcc64def0d6half_ky_plane

Last Change
Raw File

geometry_mod.F90
View Options

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
implicit none
public
contains
subroutine eval_magnetic_geometry
! evalute metrix, elements, jacobian and gradient
implicit none
!
IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
! circular model
call eval_circular_geometry
!
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
!
subroutine eval_circular_geometry
! evaluate circular geometry model
! Ref: Lapilonne et al., PoP, 2009
implicit none
REAL(dp) :: z, kx, ky
zloop: DO iz = izs,ize
z = zarray(iz)
! Metric elements
gxx(iz) = 1._dp
gxy(iz) = shear*z -eps*SIN(z)
gyy(iz) = 1._dp +(shear*z)**2 -2._dp*eps*COS(z) -2._dp*shear*eps*z*SIN(z)
! Relative strengh of radius
hatR(iz) = 1._dp + eps*cos(z)
hatB(iz) = SQRT(gxx(iz)*gyy(iz) - gxy(iz)**2)
! Jacobian
Jacobian(iz) = q0*hatR(iz)**2
! Derivative of the magnetic field strenght
gradxB(iz) = -(COS(z) + eps*SIN(z)**2)/hatB(iz)/hatB(iz)
gradzB(iz) = eps*SIN(z)
! coefficient in the front of parallel derivative
gradz_coeff(iz) = 1._dp / Jacobian(iz) / hatB(iz)
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(ikx,iky,iz) = -SIN(z)*kx -(COS(z)+shear*z*SIN(z))*ky
ENDDO
ENDDO
IF (Nky .EQ. 1) THEN ! linear 1D run we switch kx and ky for parallel opt
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
Ckxky(ikx,iky,iz) = -SIN(z)*ky -(COS(z)+shear*z*SIN(z))*kx
ENDDO
ENDDO
ENDIF
ENDDO zloop
END SUBROUTINE eval_circular_geometry
!
!--------------------------------------------------------------------------------
! UNUSED
! subroutine eval_salpha_geometry
! ! evaluate s-alpha geometry model
! implicit none
! REAL(dp) :: z, kx, ky
!
! zloop: DO iz = izs,ize
! z = zarray(iz)
! gxx(iz) = 1._dp
! gxy(iz) = shear*z
! gyy(iz) = 1._dp + (shear*z)**2
! gyz(iz) = 1._dp/eps
! gxz(iz) = 0._dp
!
! ! Relative strengh of radius
! hatR(iz) = 1._dp + eps*cos(z)
!
! ! Jacobian
! Jacobian(iz) = q0*hatR(iz)
!
! ! Relative strengh of modulus of B
! hatB(iz) = 1._dp / hatR( iz)
!
! ! Derivative of the magnetic field strenght
! gradxB(iz) = - cos( z) / hatR(iz)
! gradzB(iz) = eps * sin(z)
!
! ! Gemoetrical coefficients for the curvature operator
! ! Note: Gamma2 and Gamma3 are obtained directly form Gamma1 in the expression of the curvature operator implemented here
! !
! Gamma1(iz) = gxy(iz) * gxy(iz) - gxx(iz) * gyy(iz)
! Gamma2(iz) = gxz(iz) * gxy(iz) - gxx(iz) * gyz(iz)
! Gamma3(iz) = gxz(iz) * gyy(iz) - gxy(iz) * gyz(iz)
!
! ! Curvature operator
! DO iky = ikys, ikye
! ky = kyarray(iky)
! DO ikx= ikxl, ikxr
! kx = kxarray(ikx,iky)
! Cxy(ikx, iky, iz) = (-sin(z)*kx - (cos(z) + shear* z* sin(z))*ky) * hatB(iz) ! .. 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) = 1._dp / Jacobian(iz) / hatB(iz)
! ENDDO
!
! ! 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 iky = ikys, ikye
! ky = kyarray(iky)
! DO ikx = ikxl, ikxr
! kx = kxarray(ikx,iky)
! kperp_array(ikx, iky, iz) = sqrt( gxx(iz)*kx**2 + 2._dp* gxy(iz) * kx*ky + gyy(iz)* ky**2) !! / hatB( iz) ! there is a factor 1/B from the normalization; important to match GENE
! ENDDO
! ENDDO
! ENDDO zloop
! !
! IF( me .eq. 0 ) PRINT*, 'max kperp = ', maxval( kperp_array)
!
! END SUBROUTINE eval_salpha_geometry
!
!--------------------------------------------------------------------------------
!
end module geometry
c4science · Help