From 757a1a4ae5b671e4f117a212221061c4676a95ef Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine-hoffmann@hotmail.com>
Date: Tue, 19 Oct 2021 14:01:21 +0200
Subject: [PATCH] Correction of 3D geometry
---
src/geometry_mod.F90 | 74 +++++++++++++++++++++++++++++++++++++++---
src/moments_eq_rhs.F90 | 4 +--
src/numerics_mod.F90 | 8 ++---
3 files changed, 76 insertions(+), 10 deletions(-)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 4703bf1c..58b4e53c 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -18,15 +18,15 @@ contains
! evalute metrix, elements, jacobian and gradient
implicit none
!
- IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha geometry'
+ IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
! circular model
- call eval_s_alpha_geometry
+ call eval_circular_geometry
!
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
!
- subroutine eval_s_alpha_geometry
+ subroutine eval_circular_geometry
! evaluate circular geometry model
! Ref: Lapilonne et al., PoP, 2009
@@ -73,6 +73,72 @@ contains
ENDDO
ENDIF
ENDDO zloop
- END SUBROUTINE eval_s_alpha_geometry
+ 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
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs.F90
index ce74ae12..60699448 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs.F90
@@ -55,7 +55,7 @@ SUBROUTINE moments_eq_rhs_e
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
!! Compute moments mixing terms
! Perpendicular dynamic
@@ -199,7 +199,7 @@ SUBROUTINE moments_eq_rhs_i
ky = kyarray(iky) ! toroidal wavevector
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
!! Compute moments mixing terms
! Perpendicular dynamic
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index e8f1b486..fbbaf93a 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -82,7 +82,7 @@ SUBROUTINE evaluate_kernels
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
be_2 = kperp2 * sigmae2_taue_o2
kernel_e(ij, ikx, iky,iz) = be_2**j_int * exp(-be_2)/factj
ENDDO
@@ -95,7 +95,7 @@ SUBROUTINE evaluate_kernels
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
be_2 = kperp2 * sigmae2_taue_o2
! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
kernel_e(ijeg_e,ikx,iky,iz) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky,iz)
@@ -127,7 +127,7 @@ SUBROUTINE evaluate_kernels
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
bi_2 = kperp2 * sigmai2_taui_o2
kernel_i(ij, ikx, iky,iz) = bi_2**j_int * exp(-bi_2)/factj
ENDDO
@@ -140,7 +140,7 @@ SUBROUTINE evaluate_kernels
DO iky = ikys,ikye
ky = kyarray(iky)
DO iz = izs,ize
- kperp2= kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
+ kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
bi_2 = kperp2 * sigmai2_taui_o2
! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
kernel_i(ijeg_i,ikx,iky,iz) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky,iz)
--
GitLab