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

gyacomo/src/geometry_mod.F90
757a1a4ae5b6half_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