From 7fed0926d918453198ce8091c9e1a9f0b66b00ed Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 7 Aug 2023 17:33:13 +0200
Subject: [PATCH] Adaptation for ExB shear flow -kx is now 2D (ky dependant)
-kperp^2 array for computation reduction -routines to update the grids
---
src/collision_mod.F90 | 20 +++----
src/geometry_mod.F90 | 61 +++++++++++++---------
src/grid_mod.F90 | 104 +++++++++++++++++++++++--------------
src/inital.F90 | 18 +++----
src/model_mod.F90 | 12 +++--
src/moments_eq_rhs_mod.F90 | 8 +--
src/numerics_mod.F90 | 14 ++---
src/prec_const_mod.F90 | 2 +-
8 files changed, 139 insertions(+), 100 deletions(-)
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 5d842e4d..66df93ec 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -116,14 +116,14 @@ CONTAINS
!******************************************************************************!
SUBROUTINE compute_lenard_bernstein
USE grid, ONLY: ias,iae, ips,ipe, ijs,ije, parray, jarray, &
- ikys,ikye, ikxs,ikxe, izs,ize, kparray, ngp, ngj, ngz
+ ikys,ikye, ikxs,ikxe, izs,ize, kp2array, ngp, ngj, ngz
USE species, ONLY: sigma2_tau_o2, nu_ab
USE time_integration, ONLY: updatetlevel
USE array, ONLY: Capj
USE fields, ONLY: moments
IMPLICIT NONE
COMPLEX(xp) :: TColl_
- REAL(xp) :: j_xp, p_xp, ba_2, kp
+ REAL(xp) :: j_xp, p_xp, ba_2, kp2
INTEGER :: iz,ikx,iky,ij,ip,ia,eo,ipi,iji,izi
DO iz = izs,ize; izi = iz + ngz/2
DO ikx = ikxs, ikxe;
@@ -135,8 +135,8 @@ CONTAINS
eo = MODULO(parray(ipi),2)
p_xp = REAL(parray(ipi),xp)
j_xp = REAL(jarray(iji),xp)
- kp = kparray(iky,ikx,izi,eo)
- ba_2 = kp**2 * sigma2_tau_o2(ia) ! this is (ba/2)^2
+ kp2 = kp2array(iky,ikx,izi,eo)
+ ba_2 = kp2 * sigma2_tau_o2(ia) ! this is (ba/2)^2
!** Assembling collison operator **
! -nuee (p + 2j) Nepj
@@ -160,7 +160,7 @@ CONTAINS
USE grid, ONLY: local_na, local_np, local_nj, parray, jarray, ngp, ngj, &
ip0, ip2, ij0, ij1, ieven, &
local_nky, local_nkx, local_nz, ngz,&
- kparray
+ kp2array
USE species, ONLY: nu_ab, tau, q_tau
USE time_integration, ONLY: updatetlevel
USE array, ONLY: Capj
@@ -169,12 +169,12 @@ CONTAINS
IMPLICIT NONE
COMPLEX(xp) :: Tmp
INTEGER :: iz,ikx,iky,ij,ip,ia, ipi,iji,izi
- REAL(xp) :: j_xp, p_xp, kp_xp
+ REAL(xp) :: j_xp, p_xp, kp2_xp
DO iz = 1,local_nz
izi = iz + ngz/2
DO ikx = 1,local_nkx
DO iky = 1,local_nky
- kp_xp = kparray(iky,ikx,izi,ieven)
+ kp2_xp = kp2array(iky,ikx,izi,ieven)
DO ij = 1,local_nj
iji = ij + ngj/2
DO ip = 1,local_np
@@ -185,17 +185,17 @@ CONTAINS
j_xp = REAL(jarray(iji),xp)
!** Assembling collison operator **
IF( (p_xp .EQ. 0._xp) .AND. (j_xp .EQ. 0._xp)) THEN !Ca00
- Tmp = tau(ia)**2 * kp_xp**4*(&
+ Tmp = tau(ia)**2 * kp2_xp**2*(&
67_xp/120_xp *moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
+67_xp*SQRT2/240_xp*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel)&
-67_xp/240_xp *moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel)&
-3_xp/10_xp *q_tau(ia)*phi(iky,ikx,izi))
ELSEIF( (p_xp .EQ. 2._xp) .AND. (j_xp .EQ. 0._xp)) THEN ! Ca20
- Tmp = tau(ia) * kp_xp**2*(&
+ Tmp = tau(ia) * kp2_xp*(&
-SQRT2*twothird*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
-2._xp*twothird*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel))
ELSEIF( (p_xp .EQ. 0._xp) .AND. (j_xp .EQ. 1._xp)) THEN ! Ca01
- Tmp = tau(ia) * kp_xp**2*(&
+ Tmp = tau(ia) * kp2_xp*(&
2._xp*twothird*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
-2._xp*twothird*moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel))
ELSE
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 4d18d1c9..b20224bc 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -51,7 +51,7 @@ implicit none
! derivatives of magnetic field strength
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dBdx, dBdy, dBdz, dlnBdz
! Relative magnetic field strength
- REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, inv_hatB2 ! normalized bckg magnetic gradient amplitude and its inv squared
! Relative strength of major radius
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
@@ -71,7 +71,7 @@ implicit none
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
- eval_magnetic_geometry, set_ikx_zBC_map
+ eval_magnetic_geometry, set_ikx_zBC_map, evaluate_magn_curv
CONTAINS
@@ -164,15 +164,39 @@ CONTAINS
ERROR STOP '>> ERROR << geometry not recognized!!'
END SELECT
ENDIF
+ ! inv squared of the magnetic gradient bckg amplitude (for fast kperp update)
+ inv_hatB2 = 1/hatB/hatB
! Reset kx grid (to account for Cyq0_x0 factor)
CALL set_kxgrid(shear,Npol,Cyq0_x0,LINEARITY,N_HD)
!
! 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_
- CALL set_kparray(gxx,gxy,gyy,hatB)
+ CALL set_kparray(gxx,gxy,gyy,inv_hatB2)
+ ! Evaluate magnetic curvature operator Ckxky
+ CALL evaluate_magn_curv
+ ! coefficient in the front of parallel derivative
+ gradz_coeff = 1._xp /(jacobian*hatB)
+ ! d/dz(ln B) to correspond to the formulation in Hoffmann et al. 2023
+ dlnBdz = dBdz/hatB
+ !
+ ! set the mapping for parallel boundary conditions
+ CALL set_ikx_zBC_map
+ !
+ ! Compute the inverse z integrated Jacobian (useful for flux averaging)
+ integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
+ CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
+ iInt_Jacobian = 1._xp/iInt_Jacobian ! reverse it
+ END SUBROUTINE eval_magnetic_geometry
+
+ SUBROUTINE evaluate_magn_curv
+ USE grid, ONLY: local_nkx, local_nky, local_nz, ngz,&
+ kxarray, kyarray, nzgrid
+ IMPLICIT NONE
+ REAL(xp) :: kx,ky
+ real(xp) :: Cx, Cy, g_xz, g_yz, g_zz
+ INTEGER :: eo,iz,iky,ikx
DO eo = 1,nzgrid
- ! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = 1,local_nz+ngz
! !covariant metric
g_xz = jacobian(iz,eo)**2*(gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo))
@@ -186,25 +210,13 @@ CONTAINS
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx= 1,local_nkx
- kx = kxarray(ikx)
+ kx = kxarray(iky,ikx)
Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/C_xy
ENDDO
ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._xp /(jacobian(iz,eo)*hatB(iz,eo))
- ! d/dz(ln B) to correspond to the formulation in paper 2023
- dlnBdz(iz,eo) = dBdz(iz,eo)/hatB(iz,eo)
ENDDO
ENDDO
- !
- ! set the mapping for parallel boundary conditions
- CALL set_ikx_zBC_map
- !
- ! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
- CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
- iInt_Jacobian = 1._xp/iInt_Jacobian ! reverse it
- END SUBROUTINE eval_magnetic_geometry
+ END SUBROUTINE
!
!--------------------------------------------------------------------------------
!
@@ -382,10 +394,10 @@ CONTAINS
! Check if it has to be connected to the otherside of the kx grid
if(ikx_zBC_L(iky,ikx) .LE. 0) ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + total_nkx
! Check if it points out of the kx domain
- ! print*, kxarray(ikx), shift, kx_min
- IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain
- ! print*,( ((kxarray(ikx) - shift)-kx_min) .LT. EPSILON(kx_min) )
- ! print*, kxarray(ikx)-kx_min, EPSILON(kx_min)
+ ! print*, kxarray(iky,ikx), shift, kx_min
+ IF( (kxarray(iky,ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain
+ ! print*,( ((kxarray(iky,ikx) - shift)-kx_min) .LT. EPSILON(kx_min) )
+ ! print*, kxarray(iky,ikx)-kx_min, EPSILON(kx_min)
! IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN
SELECT CASE(parallel_bc)
CASE ('dirichlet','disconnected') ! connected to 0
@@ -420,7 +432,7 @@ CONTAINS
! Check if it has to be connected to the otherside of the kx grid
if(ikx_zBC_R(iky,ikx) .GT. total_nkx) ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - total_nkx
! Check if it points out of the kx domain
- IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
+ IF( (kxarray(iky,ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
! IF( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain
SELECT CASE(parallel_bc)
CASE ('dirichlet','disconnected') ! connected to 0
@@ -428,7 +440,7 @@ CONTAINS
CASE ('periodic') ! connected to itself as for shearless
ikx_zBC_R(iky,ikx) = ikx
CASE ('cyclic')
- ! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max
+ ! write(*,*) 'check',ikx,iky, kxarray(iky,ikx) + shift, '>', kx_max
ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,total_nkx)+1
CASE DEFAULT
ERROR STOP "The parallel BC is not recognized (dirichlet, periodic, cyclic or disconnected)"
@@ -517,6 +529,7 @@ END SUBROUTINE set_ikx_zBC_map
ALLOCATE( dBdy(local_nz+ngz,nzgrid))
ALLOCATE( dBdz(local_nz+ngz,nzgrid))
ALLOCATE( dlnBdz(local_nz+ngz,nzgrid))
+ ALLOCATE( inv_hatB2(local_nz+ngz,nzgrid))
ALLOCATE( hatB(local_nz+ngz,nzgrid))
ALLOCATE( hatR(local_nz+ngz,nzgrid))
ALLOCATE( hatZ(local_nz+ngz,nzgrid))
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index dbeda1d1..5e52854b 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -25,11 +25,13 @@ MODULE grid
! Grid arrays
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: parray, parray_full
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: jarray, jarray_full
- REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray, kxarray_full
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray ! ExB shear makes it ky dependant
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray_full
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full
REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full
REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp
+ REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kp2array !kperp^2
! Kronecker delta for p=0, p=1, p=2, j=0, j=1
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kroneck_p0, kroneck_p1, kroneck_p2, kroneck_p3
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kroneck_j0, kroneck_j1
@@ -112,6 +114,7 @@ MODULE grid
PUBLIC :: set_grids, set_kxgrid, set_kparray
PUBLIC :: grid_readinputs, grid_outputinputs
PUBLIC :: bar
+ PUBLIC :: update_grids ! Update the kx and kperp grids for ExB shear flow
! Precomputations
real(xp), PUBLIC, PROTECTED :: pmax_xp, jmax_xp
@@ -214,7 +217,7 @@ CONTAINS
local_nkx_ptr = ikxe - ikxs + 1
local_nkx = ikxe - ikxs + 1
local_nkx_offset = ikxs - 1
- ALLOCATE(kxarray(local_nkx))
+ ALLOCATE(kxarray(local_nky,local_Nkx))
ALLOCATE(kxarray_full(total_nkx))
ALLOCATE(AA_x(local_nkx))
!!---------------- PARALLEL Z GRID (parallelized)
@@ -263,6 +266,7 @@ CONTAINS
ENDDO
!!---------------- Kperp grid
CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
+ CALL allocate_array(kp2array, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
END SUBROUTINE
!!!! GRID REPRESENTATION
! We define the grids that contain ghosts (p,j,z) with indexing from 1 to Nlocal + nghost
@@ -474,7 +478,7 @@ CONTAINS
REAL(xp), INTENT(IN) :: shear, Npol, Cyq0_x0
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
- INTEGER :: ikx
+ INTEGER :: ikx, iky
REAL(xp):: Lx_adapted
IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
@@ -499,18 +503,20 @@ CONTAINS
kx_min = MINVAL(kxarray_full)!-kx_max+deltakx
! Set local grid (not parallelized so same as full one)
local_kxmax = 0._xp
- DO ikx = 1,local_nkx
- kxarray(ikx) = kxarray_full(ikx+local_nkx_offset)
- ! Finding kx=0
- IF (kxarray(ikx) .EQ. 0) THEN
- ikx0 = ikx
- contains_kx0 = .true.
- ENDIF
- ! Finding local kxmax
- IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
- local_kxmax = ABS(kxarray(ikx))
- ikx_max = ikx
- ENDIF
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ kxarray(iky,ikx) = kxarray_full(ikx+local_nkx_offset)
+ ! Finding kx=0
+ IF (kxarray(1,ikx) .EQ. 0) THEN
+ ikx0 = ikx
+ contains_kx0 = .true.
+ ENDIF
+ ! Finding local kxmax
+ IF (ABS(kxarray(iky,ikx)) .GT. local_kxmax) THEN
+ local_kxmax = ABS(kxarray(iky,ikx))
+ ikx_max = ikx
+ ENDIF
+ END DO
END DO
ELSE ! Odd number of kx (-2 -1 0 1 2)
ERROR STOP "Gyacomo is safer with an even Kx number"
@@ -518,13 +524,15 @@ CONTAINS
! Orszag 2/3 filter
two_third_kxmax = 2._xp/3._xp*kx_max;
! Antialiasing filter
- DO ikx = 1,local_nkx
- IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. &
- (kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_x(ikx) = 1._xp;
- ELSE
- AA_x(ikx) = 0._xp;
- ENDIF
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ IF ( ((kxarray(iky,ikx) .GT. -two_third_kxmax) .AND. &
+ (kxarray(iky,ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
+ AA_x(ikx) = 1._xp;
+ ELSE
+ AA_x(ikx) = 0._xp;
+ ENDIF
+ END DO
END DO
! For hyperdiffusion
IF(LINEARITY.EQ.'linear') THEN
@@ -572,19 +580,6 @@ CONTAINS
! Find local extrema
local_zmax(eo) = zarray(local_nz+ngz/2,eo)
local_zmin(eo) = zarray(1+ngz/2,eo)
- ! Periodic z \in (-pi pi-dz)
- ! DO ig = 1,ngz/2 ! first ghost cells
- ! iglob = ig+local_nz_offset-ngz/2
- ! IF (iglob .LE. 0) &
- ! iglob = iglob + total_nz
- ! zarray(ig,eo) = zarray_full(iglob)
- ! ENDDO
- ! DO ig = local_nz+ngz/2,local_nz+ngz ! last ghost cells
- ! iglob = ig+local_nz_offset-ngz/2
- ! IF (iglob .GT. total_nz) &
- ! iglob = iglob - total_nz
- ! zarray(ig,eo) = zarray_full(iglob)
- ! ENDDO
! continue in z<pi and z>pi
DO ig = 1,ngz/2
zarray(local_nz+ngz/2+ig,eo) = zarray(local_nz+ngz/2,eo) + ig*deltaz
@@ -611,9 +606,9 @@ CONTAINS
ENDIF
END SUBROUTINE set_zgrid
- SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
+ SUBROUTINE set_kparray(gxx, gxy, gyy, inv_hatB2)
IMPLICIT NONE
- REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,hatB
+ REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
INTEGER :: eo,iz,iky,ikx
REAL(xp) :: kx, ky
DO eo = 1,nzgrid
@@ -621,11 +616,12 @@ CONTAINS
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx = 1,local_nkx
- kx = kxarray(ikx)
+ kx = kxarray(iky,ikx)
! there is a factor 1/B from the normalization; important to match GENE
! this factor comes from $b_a$ argument in the Bessel. Kperp is not used otherwise.
- kparray(iky, ikx, iz, eo) = &
- SQRT( gxx(iz,eo)*kx**2 + 2._xp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/ hatB(iz,eo)
+ kp2array(iky, ikx, iz, eo) = &
+ (gxx(iz,eo)*kx**2 + 2._xp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)*inv_hatB2(iz,eo)
+ kparray(iky, ikx, iz, eo) = SQRT(kp2array(iky, ikx, iz, eo))
ENDDO
ENDDO
ENDDO
@@ -633,6 +629,34 @@ CONTAINS
two_third_kpmax = 2._xp/3._xp * MAXVAL(kparray)
END SUBROUTINE
+ SUBROUTINE update_grids (sky,gxx,gxy,gyy,inv_hatB2)
+ IMPLICIT NONE
+ REAL(xp), DIMENSION(local_nky),INTENT(IN) :: sky ! ExB grid shift
+ REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
+ INTEGER :: eo,iz,iky,ikx
+ REAL(xp) :: kx, ky, skp2
+ ! Update the kx grid
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ kxarray(iky,ikx) = kxarray(iky,ikx) + sky(iky)
+ ENDDO
+ ENDDO
+ ! Update the kperp grid
+ DO eo = 1,nzgrid
+ DO iz = 1,local_nz+ngz
+ DO iky = 1,local_nky
+ ky = kyarray(iky)
+ DO ikx = 1,local_nkx
+ kx = kxarray(iky,ikx)
+ ! shift in kp obtained from kp2* = kp2 + skp2
+ skp2 = sky(iky)*inv_hatB2(iz,eo)*(gxx(iz,eo)*(2._xp*kx+sky(iky)) +gxy(iz,eo)*ky)
+ kp2array(iky, ikx, iz, eo) = kp2array(iky, ikx, iz, eo) + skp2
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ END SUBROUTINE
+
SUBROUTINE grid_outputinputs(fid)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach, creatd
diff --git a/src/inital.F90 b/src/inital.F90
index cfc12ecb..8871714a 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -308,10 +308,10 @@ SUBROUTINE init_phi_ppj
INTEGER :: ikx, iky, iz
amp = 1.0_xp
!**** ppj initialization *******************************************
- DO ikx=1,total_nkx
- kx = kxarray(ikx)
- DO iky=1,local_nky
- ky = kyarray(iky)
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
DO iz=1,local_nz+ngz
z = zarray(iz,ieven)
IF (ky .NE. 0) THEN
@@ -381,7 +381,7 @@ SUBROUTINE initialize_blob
DO iz=1+ngz/2,local_nz+ngz/2
z = zarray(iz,ieven)
DO ikx=1,total_nkx
- kx = kxarray(ikx) + z*shear*ky
+ kx = kxarray(iky,ikx) + z*shear*ky
DO ip=1+ngp/2,local_np+ngp/2
p = parray(ip)
DO ij=1+ngj/2,local_nj+ngj/2
@@ -437,10 +437,10 @@ SUBROUTINE init_ppj
DO ip=1,local_np+ngp
DO ij=1,local_nj+ngj
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
- DO ikx=1,total_nkx
- kx = kxarray(ikx)
- DO iky=1,local_nky
- ky = kyarray(iky)
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
DO iz=1,local_nz+ngz
z = zarray(iz,ieven)
IF (kx .EQ. 0) THEN
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 4f4974ca..abd264ee 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -17,14 +17,15 @@ MODULE model
CHARACTER(len=16), &
PUBLIC, PROTECTED :: HYP_V = 'hypcoll' ! hyperdiffusion model for velocity space ('none','hypcoll','dvpar4')
INTEGER, PUBLIC, PROTECTED :: Na = 1 ! number of evolved species
+ LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
+ LOGICAL, PUBLIC :: ADIAB_I = .false. ! adiabatic ion model
+ REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
REAL(xp), PUBLIC, PROTECTED :: nu = 0._xp ! collision frequency parameter
REAL(xp), PUBLIC, PROTECTED :: k_gB = 1._xp ! Magnetic gradient strength (L_ref/L_gB)
REAL(xp), PUBLIC, PROTECTED :: k_cB = 1._xp ! Magnetic curvature strength (L_ref/L_cB)
REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
- LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
- LOGICAL, PUBLIC :: ADIAB_I = .false. ! adiabatic ion model
- REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
+ REAL(xp), PUBLIC, PROTECTED :: ExBrate = 0._xp ! ExB background shearing rate (radially constant shear flow)
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift
@@ -46,8 +47,9 @@ CONTAINS
IMPLICIT NONE
NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
- mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, Na,&
- nu, k_gB, k_cB, lambdaD, beta, ADIAB_E, ADIAB_I, tau_i
+ Na, ADIAB_E, ADIAB_I, tau_i, &
+ mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, &
+ nu, k_gB, k_cB, lambdaD, beta, ExBrate
READ(lu_in,model_par)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index ea1691cb..f7609cb1 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -10,7 +10,7 @@ CONTAINS
USE grid, ONLY: local_na, local_np, local_nj, local_nkx, local_nky, local_nz,&
nzgrid,pp2,ngp,ngj,ngz,&
diff_dz_coeff,diff_kx_coeff,diff_ky_coeff,diff_p_coeff,diff_j_coeff,&
- parray,jarray,kxarray, kyarray, kparray
+ parray,jarray,kxarray, kyarray, kp2array
USE basic
USE closure, ONLY: evolve_mom
USE prec_const
@@ -35,10 +35,10 @@ CONTAINS
z:DO iz = 1,local_nz
izi = iz + ngz/2
x:DO ikx = 1,local_nkx
- kx = kxarray(ikx) ! radial wavevector
- i_kx = imagu * kx ! radial derivative
y:DO iky = 1,local_nky
+ kx = kxarray(iky,ikx) ! radial wavevector
ky = kyarray(iky) ! binormal wavevector
+ i_kx = imagu * kx ! radial derivative
i_ky = imagu * ky ! binormal derivative
phikykxz = phi(iky,ikx,izi)
psikykxz = psi(iky,ikx,izi)
@@ -50,7 +50,7 @@ CONTAINS
ipi = ip+ngp/2
p_int = parray(ipi) ! Hermite degree
eo = min(nzgrid,MODULO(p_int,2)+1) ! Indicates if we are on odd or even z grid
- kperp2= kparray(iky,ikx,izi,eo)**2
+ kperp2= kp2array(iky,ikx,izi,eo)
! Species loop
a:DO ia = 1,local_na
Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 0224a515..6ec4d1de 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -74,7 +74,7 @@ SUBROUTINE evaluate_kernels
USE basic
USE prec_const, ONLY: xp
USE array, ONLY : kernel!, HF_phi_correction_operator
- USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kparray,&
+ USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kp2array,&
nzgrid
USE species, ONLY : sigma2_tau_o2
USE model, ONLY : ADIAB_I, tau_i
@@ -91,7 +91,7 @@ DO ia = 1,local_na
DO ij = 1,local_nj + ngj
j_int = jarray(ij)
j_xp = REAL(j_int,xp)
- y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
+ y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo)
IF(j_int .LT. 0) THEN !ghosts values
kernel(ia,ij,iky,ikx,iz,eo) = 0._xp
ELSE
@@ -114,7 +114,7 @@ DO ia = 1,local_na
DO ij = 1,local_nj + ngj
j_int = jarray(ij)
j_xp = REAL(j_int,xp)
- y_ = sigma_i**2*tau_i/2._xp * kparray(iky,ikx,iz,eo)**2
+ y_ = sigma_i**2*tau_i/2._xp * kp2array(iky,ikx,iz,eo)
IF(j_int .LT. 0) THEN !ghosts values
kernel_i(ij,iky,ikx,iz,eo) = 0._xp
ELSE
@@ -174,7 +174,7 @@ SUBROUTINE evaluate_poisson_op
kxloop: DO ikx = 1,local_nkx
kyloop: DO iky = 1,local_nky
zloop: DO iz = 1,local_nz
- IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ IF( (kxarray(iky,ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
inv_poisson_op(iky, ikx, iz) = 0._xp
inv_pol_ion (iky, ikx, iz) = 0._xp
ELSE
@@ -217,7 +217,7 @@ SUBROUTINE evaluate_ampere_op
USE prec_const, ONLY : xp
USE array, ONLY : kernel, inv_ampere_op
USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, ngz, total_nj, ngj,&
- kparray, kxarray, kyarray, SOLVE_AMPERE, iodd
+ kp2array, kxarray, kyarray, SOLVE_AMPERE, iodd
USE model, ONLY : beta, ADIAB_I
USE species, ONLY : q, sigma
USE geometry, ONLY : hatB
@@ -232,8 +232,8 @@ SUBROUTINE evaluate_ampere_op
x:DO ikx = 1,local_nkx
y:DO iky = 1,local_nky
z:DO iz = 1,local_nz
- kperp2 = kparray(iky,ikx,iz+ngz/2,iodd)**2
- IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ kperp2 = kp2array(iky,ikx,iz+ngz/2,iodd)
+ IF( (kxarray(iky,ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
inv_ampere_op(iky, ikx, iz) = 0._xp
ELSE
sum_jpol = 0._xp
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index ce1efa47..f8fedacd 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -78,7 +78,7 @@ MODULE prec_const
! dp_r = range(b)
! dp_p = precision(b)
- REAL(xp) :: c = 1_xp
+ REAL(xp) :: c = 1._xp
xp_r = range(c)
xp_p = precision(c)
--
GitLab