diff --git a/src/array_mod.F90 b/src/array_mod.F90
index 1cb6d141da980927dcce532bea67e626bbc25672..16a1bdf32392a716c9321f4b8a3b4ab113f23f14 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -6,6 +6,9 @@ MODULE array
! Arrays to store the rhs, for time integration (ip,ij,ikx,iky,iz,updatetlevel)
COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_e
COMPLEX(dp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: moments_rhs_i
+ ! Non linear term array (ip,ij,ikx,iky,iz)
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron
+ COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sipj ! ion
! To load collision matrix (ip1,ij1,ip2,ij2)
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ceepj, CeipjT
@@ -26,17 +29,31 @@ MODULE array
REAL(dp), DIMENSION(:), ALLOCATABLE :: xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
- ! Curvature array
- COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ckxky
-
+ ! Geoemtrical operators
+ ! Curvature
+ REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z
+ ! Jacobian
+ REAL(dp), DIMENSION(:), ALLOCATABLE :: Jacobian ! dimensions: z
+ ! Metric
+ REAL(dp), DIMENSION(:), ALLOCATABLE :: gxx, gxy, gyy, gxz, gyz
+ ! derivatives of magnetic field strength
+ REAL(dp), DIMENSION(:), allocatable :: gradzB ! dimensions: z
+ REAL(dp), DIMENSION(:), allocatable :: gradxB
+ ! Relative magnetic field strength
+ REAL(dp), DIMENSION(:), allocatable :: hatB
+ ! Relative strength of major radius
+ REAL(dp), DIMENSION(:), allocatable :: hatR
+ ! Geometrical factors
+ REAL(dp), DIMENSION(:), allocatable :: Gamma1
+ REAL(dp), DIMENSION(:), allocatable :: Gamma2
+ REAL(dp), DIMENSION(:), allocatable :: Gamma3
+ ! Some geometrical coefficients
+ REAL(dp), DIMENSION(:) , allocatable :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Kernel function evaluation (ij,ikx,iky)
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_e
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: kernel_i
- ! Non linear term array (ip,ij,ikx,iky,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sepj ! electron
- COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: Sipj ! ion
-
+ !! Diagnostics
! Gyrocenter density for electron and ions (ikx,iky,iz)
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ne00
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Ni00
diff --git a/src/auxval.F90 b/src/auxval.F90
index f1c3fd6935474a90b282e8297dff3cc3241ab2aa..22d63d44a7281ab7f57f2aaf524528ef3cf0700f 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -8,6 +8,7 @@ subroutine auxval
USE fourier, ONLY: init_grid_distr_and_plans, alloc_local_1, alloc_local_2
use prec_const
USE numerics
+ USE geometry
IMPLICIT NONE
INTEGER :: irows,irowe, irow, icol, i_
@@ -29,7 +30,9 @@ subroutine auxval
CALL memory ! Allocate memory for global arrays
- CALL lin_coeff_and_geometry ! precompute coeff for lin equation and geometry
+ CALL eval_magnetic_geometry ! precompute coeff for lin equation
+
+ CALL compute_lin_coeff ! precompute coeff for lin equation and geometry
CALL evaluate_kernels ! precompute the kernels
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..2ebda8f58f3260142c1a780927c2abe9ef7f2d22
--- /dev/null
+++ b/src/geometry_mod.F90
@@ -0,0 +1,95 @@
+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) :: shear = 0._dp
+ REAL(dp) :: epsilon = 0.18_dp
+
+ ! Metric elements
+
+ DO iz = izs,ize
+ gxx(iz) = 1._dp
+ gxy(iz) = shear*zarray(iz) - epsilon*sin(zarray(iz))
+ gyy(iz) = 1._dp + (shear*zarray(iz))**2 &
+ - 2._dp * epsilon *COS(zarray(iz)) - 2._dp*shear * epsilon * zarray(iz)*SIN(zarray(iz))
+ gxz( iz) = - SIN(zarray(iz))
+ gyz( iz) = ( 1._dp - 2._dp * epsilon *COS(zarray( iz)) - epsilon*shear * zarray( iz) * SIN(zarray(iz)) ) /epsilon
+ ENDDO
+
+ ! Relative strengh of radius
+ DO iz = izs,ize
+ hatR(iz) = 1._dp + epsilon*cos(zarray(iz))
+ ENDDO
+
+ DO iz = izs, ize
+ hatB( iz) = sqrt(gxx(iz) * gyy(iz) - gxy(iz)* gxy(iz))
+ ENDDO
+
+
+ ! Jacobian
+ DO iz = izs,ize
+ Jacobian(iz) = q0*hatR(iz)*hatR(iz)
+ ENDDO
+
+ ! Derivative of the magnetic field strenght
+ DO iz = izs, ize
+ gradxB(iz) = - ( COS( zarray(iz)) + epsilon* SIN( zarray(iz)) * SIN(zarray(iz) )) / hatB(iz)/hatB(iz)
+ gradzB( iz) = epsilon * SIN(zarray(iz)) *( 1._dp - epsilon*COS(zarray(iz)) ) / hatB(iz)/hatB(iz)
+ ENDDO
+
+ ! Gemoetrical coefficients for the curvature operator
+ ! Note: Gamma2 and Gamma3 are obtained directly form Gamma1 in the expression of the curvature operator implemented here
+ DO iz = izs, ize
+ 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)
+ ENDDO
+
+ ! Curvature operator
+ DO iz = izs, ize
+ DO iky = ikys, ikye
+ DO ikx= ikxs, ikxe
+ Ckxky( ikx, iky, iz) = Gamma1(iz)/hatB(iz)*((kxarray(ikx) &
+ + shear*zarray(iz)*kyarray(iky)) * gradzB(iz)/epsilon &
+ - gradxB(iz)* kyarray(iky))
+ ENDDO
+ ENDDO
+ ENDDO
+
+ ! coefficient in the front of parallel derivative
+ DO iz = izs, ize
+ gradz_coeff( iz) = 1._dp / Jacobian( iz) / hatB( iz)
+ ENDDO
+
+ END SUBROUTINE eval_circular_geometry
+
+end module geometry
diff --git a/src/lin_coeff_and_geometry.F90 b/src/lin_coeff_and_geometry.F90
deleted file mode 100644
index f2736623f03389363fe70c93f897242bba01a66d..0000000000000000000000000000000000000000
--- a/src/lin_coeff_and_geometry.F90
+++ /dev/null
@@ -1,99 +0,0 @@
-SUBROUTINE lin_coeff_and_geometry
- USE array, ONLY: xnepj, xnepp1j, xnepm1j, xnepp2j, xnepm2j, xnepjp1, xnepjm1, &
- xnipj, xnipp1j, xnipm1j, xnipp2j, xnipm2j, xnipjp1, xnipjm1, &
- xphij, xphijp1, xphijm1, Ckxky
- USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, eta_T, eta_n
- USE prec_const
- USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
- ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i,&
- kxarray, kyarray, zarray, &
- ikx,iky,iz, ikxs,ikxe, ikys,ikye, izs,ize
- IMPLICIT NONE
- INTEGER :: p_int, j_int ! polynom. degrees
- REAL(dp) :: p_dp, j_dp
- REAL(dp) :: kx, ky, z
-
- !! Electrons linear coefficients for moment RHS !!!!!!!!!!
- DO ip = ips_e, ipe_e
- p_int= parray_e(ip) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- DO ij = ijs_e, ije_e
- j_int= jarray_e(ij) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnepj(ip,ij) = taue_qe * 2._dp * (p_dp + j_dp + 1._dp)
- ENDDO
- ENDDO
- DO ip = ips_e, ipe_e
- p_int= parray_e(ip) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- xnepp1j(ip) = sqrtTaue_qe * SQRT(p_dp + 1_dp)
- xnepm1j(ip) = sqrtTaue_qe * SQRT(p_dp)
- xnepp2j(ip) = taue_qe * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
- xnepm2j(ip) = taue_qe * SQRT(p_dp * (p_dp - 1._dp))
- ENDDO
- DO ij = ijs_e, ije_e
- j_int= jarray_e(ij) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnepjp1(ij) = -taue_qe * (j_dp + 1._dp)
- xnepjm1(ij) = -taue_qe * j_dp
- ENDDO
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !! Ions linear coefficients for moment RHS !!!!!!!!!!
- DO ip = ips_i, ipe_i
- p_int= parray_i(ip) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- DO ij = ijs_i, ije_i
- j_int= jarray_i(ij) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnipj(ip,ij) = taui_qi * 2._dp * (p_dp + j_dp + 1._dp)
- ENDDO
- ENDDO
- DO ip = ips_i, ipe_i
- p_int= parray_i(ip) ! Hermite degree
- p_dp = REAL(p_int,dp) ! REAL of Hermite degree
- xnipp1j(ip) = sqrtTaui_qi * SQRT(p_dp + 1._dp)
- xnipm1j(ip) = sqrtTaui_qi * SQRT(p_dp)
- xnipp2j(ip) = taui_qi * SQRT((p_dp + 1._dp) * (p_dp + 2._dp))
- xnipm2j(ip) = taui_qi * SQRT(p_dp * (p_dp - 1._dp))
- ENDDO
- DO ij = ijs_i, ije_i
- j_int= jarray_i(ij) ! Laguerre degree
- j_dp = REAL(j_int,dp) ! REAL of Laguerre degree
- xnipjp1(ij) = -taui_qi * (j_dp + 1._dp)
- xnipjm1(ij) = -taui_qi * j_dp
- ENDDO
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !! ES linear coefficients for moment RHS !!!!!!!!!!
- DO ip = ips_i, ipe_i
- p_int= parray_i(ip) ! Hermite degree
- DO ij = ijs_i, ije_i
- j_int= jarray_i(ij) ! REALof Laguerre degree
- j_dp = REAL(j_int,dp) ! REALof Laguerre degree
- !! Electrostatic potential pj terms
- IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij(ip,ij) = eta_n + 2.*j_dp*eta_T
- xphijp1(ip,ij) =-eta_T*(j_dp+1._dp)
- xphijm1(ip,ij) =-eta_T* j_dp
- ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij(ip,ij) = eta_T/SQRT2
- xphijp1(ip,ij) = 0._dp; xphijm1(ip,ij) = 0._dp;
- ELSE
- xphij(ip,ij) = 0._dp; xphijp1(ip,ij) = 0._dp
- xphijm1(ip,ij) = 0._dp;
- ENDIF
- ENDDO
- ENDDO
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !! Curvature and geometric coefficients !!!!!!!!!!
- DO iz = izs,ize
- z = zarray(iz)
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx) ! Poloidal wavevector
- DO iky = ikys,ikye
- ky = kyarray(iky) ! Toroidal wavevector
- Ckxky(ikx,iky,iz) = SIN(z)*imagu*kx +COS(z)*imagu*ky
- ENDDO
- ENDDO
- ENDDO
-
-END SUBROUTINE lin_coeff_and_geometry
diff --git a/src/memory.F90 b/src/memory.F90
index 6c5f3fe1051209ef813f9b4d082a02426fa2d4dc..dab37da4e1edd0056f0d168d767a81364e016761 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -73,6 +73,20 @@ SUBROUTINE memory
! Curvature and geometry
CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(Jacobian,izs,ize)
+ CALL allocate_array(gxx, izs,ize)
+ CALL allocate_array(gxy, izs,ize)
+ CALL allocate_array(gyy, izs,ize)
+ CALL allocate_array(gyz, izs,ize)
+ CALL allocate_array(gxz, izs,ize)
+ CALL allocate_array(gradzB,izs,ize)
+ CALL allocate_array(gradxB,izs,ize)
+ CALL allocate_array(hatB,izs,ize)
+ CALL allocate_array(hatR,izs,ize)
+ CALL allocate_array(Gamma1, izs,ize)
+ call allocate_array(Gamma2, izs,ize)
+ call allocate_array(Gamma3, izs, ize)
+ call allocate_array(gradz_coeff, izs, ize)
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
diff --git a/src/numerics.F90 b/src/numerics.F90
deleted file mode 100644
index b7870565722df3157fb0148d9c086924dff076ac..0000000000000000000000000000000000000000
--- a/src/numerics.F90
+++ /dev/null
@@ -1,150 +0,0 @@
-MODULE numerics
- USE basic
- USE prec_const
- USE grid
- USE utility
- USE coeff
- implicit none
-
- PUBLIC :: compute_derivatives, build_dnjs_table, evaluate_kernels
-
-CONTAINS
-
-! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
-SUBROUTINE compute_derivatives
-
-END SUBROUTINE compute_derivatives
-
-!******************************************************************************!
-!!!!!!! Build the Laguerre-Laguerre coupling coefficient table for nonlin
-!******************************************************************************!
-SUBROUTINE build_dnjs_table
- USE array, Only : dnjs
- USE coeff
- IMPLICIT NONE
-
- INTEGER :: in, ij, is, J
- INTEGER :: n_, j_, s_
-
- J = max(jmaxe,jmaxi)
-
- DO in = 1,J+1 ! Nested dependent loops to make benefit from dnjs symmetry
- n_ = in - 1
- DO ij = in,J+1
- j_ = ij - 1
- DO is = ij,J+1
- s_ = is - 1
-
- dnjs(in,ij,is) = TO_DP(ALL2L(n_,j_,s_,0))
- ! By symmetry
- dnjs(in,is,ij) = dnjs(in,ij,is)
- dnjs(ij,in,is) = dnjs(in,ij,is)
- dnjs(ij,is,in) = dnjs(in,ij,is)
- dnjs(is,ij,in) = dnjs(in,ij,is)
- dnjs(is,in,ij) = dnjs(in,ij,is)
- ENDDO
- ENDDO
- ENDDO
-END SUBROUTINE build_dnjs_table
-!******************************************************************************!
-
-!******************************************************************************!
-!!!!!!! Evaluate the kernels once for all
-!******************************************************************************!
-SUBROUTINE evaluate_kernels
- USE basic
- USE array, Only : kernel_e, kernel_i
- USE grid
- use model, ONLY : tau_e, tau_i, sigma_e, sigma_i, q_e, q_i, lambdaD, CLOS, sigmae2_taue_o2, sigmai2_taui_o2
- IMPLICIT NONE
-
- REAL(dp) :: factj, j_dp, j_int
- REAL(dp) :: be_2, bi_2, alphaD
- REAL(dp) :: kx, ky, kperp2
-
- !!!!! Electron kernels !!!!!
- !Precompute species dependant factors
- factj = 1.0 ! Start of the recursive factorial
- DO ij = 1, jmaxe+1
- j_int = jarray_e(ij)
- j_dp = REAL(j_int,dp) ! REAL of degree
-
- ! Recursive factorial
- IF (j_dp .GT. 0) THEN
- factj = factj * j_dp
- ELSE
- factj = 1._dp
- ENDIF
-
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
-
- be_2 = (kx**2 + ky**2) * sigmae2_taue_o2
-
- kernel_e(ij, ikx, iky) = be_2**j_int * exp(-be_2)/factj
-
- ENDDO
- ENDDO
- ENDDO
- ! Kernels closure
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- be_2 = (kx**2 + ky**2) * sigmae2_taue_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj (/!\ ij = j+1)
- kernel_e(ijeg_e,ikx,iky) = be_2/(real(ijeg_e-1,dp))*kernel_e(ije_e,ikx,iky)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( be_2 .NE. 0 ) THEN
- kernel_e(ijsg_e,ikx,iky) = (real(ijsg_e,dp))/be_2*kernel_e(ijs_e,ikx,iky)
- ELSE
- kernel_e(ijsg_e,ikx,iky) = 0._dp
- ENDIF
- ENDDO
- ENDDO
-
- !!!!! Ion kernels !!!!!
- factj = 1.0 ! Start of the recursive factorial
- DO ij = 1, jmaxi+1
- j_int = jarray_e(ij)
- j_dp = REAL(j_int,dp) ! REAL of degree
-
- ! Recursive factorial
- IF (j_dp .GT. 0) THEN
- factj = factj * j_dp
- ELSE
- factj = 1._dp
- ENDIF
-
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2
- kernel_i(ij, ikx, iky) = bi_2**j_int * exp(-bi_2)/factj
-
- ENDDO
- ENDDO
- ENDDO
- ! Kernels closure
- DO ikx = ikxs,ikxe
- kx = kxarray(ikx)
- DO iky = ikys,ikye
- ky = kyarray(iky)
- bi_2 = (kx**2 + ky**2) * sigmai2_taui_o2
- ! Kernel ghost + 1 with Kj+1 = y/(j+1) Kj
- kernel_i(ijeg_i,ikx,iky) = bi_2/(real(ijeg_i-1,dp))*kernel_i(ije_i,ikx,iky)
- ! Kernel ghost - 1 with Kj-1 = j/y Kj(careful about the kperp=0)
- IF ( bi_2 .NE. 0 ) THEN
- kernel_i(ijsg_i,ikx,iky) = (real(ijsg_i,dp))/bi_2*kernel_i(ijs_i,ikx,iky)
- ELSE
- kernel_i(ijsg_i,ikx,iky) = 0._dp
- ENDIF
- ENDDO
- ENDDO
-END SUBROUTINE evaluate_kernels
-!******************************************************************************!
-
-END MODULE numerics