diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index 29c08f7a5f31647569c7824739e9bce09ba0ddf4..267451912d43058f65ebd58b3202fcf5807c0344 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -216,7 +216,7 @@ END SUBROUTINE interp_z
SUBROUTINE simpson_rule_z(local_nz,zweights_SR,f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
- USE prec_const, ONLY: xp, onethird
+ USE prec_const, ONLY: xp, onethird, mpi_xp_c
USE parallel, ONLY: num_procs_z, rank_z, comm_z, manual_0D_bcast
USE mpi
implicit none
@@ -238,12 +238,12 @@ SUBROUTINE simpson_rule_z(local_nz,zweights_SR,f,intf)
IF (num_procs_z .GT. 1) THEN
!! Everyone sends its local_sum to root = 0
IF (rank_z .NE. root) THEN
- CALL MPI_SEND(buffer, 1 , MPI_DOUBLE_COMPLEX, root, 5678, comm_z, ierr)
+ CALL MPI_SEND(buffer, 1 , mpi_xp_c, root, 5678, comm_z, ierr)
ELSE
! Recieve from all the other processes
DO i_ = 0,num_procs_z-1
IF (i_ .NE. rank_z) &
- CALL MPI_RECV(buffer, 1 , MPI_DOUBLE_COMPLEX, i_, 5678, comm_z, MPI_STATUS_IGNORE, ierr)
+ CALL MPI_RECV(buffer, 1 , mpi_xp_c, i_, 5678, comm_z, MPI_STATUS_IGNORE, ierr)
intf = intf + buffer
ENDDO
ENDIF
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 9a80a8bdb57cd8bc8c2cfef5e0da1510efd3095f..7baace2eecd84016bf6d25c2d7de65e2549c4eb3 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -38,21 +38,20 @@ SUBROUTINE set_closure_model
CASE('truncation')
DO ip = 1,local_np+ngp
DO ij = 1, local_nj+ngj
- evolve_mom(ip,ij) = (parray(ip).GE.0) .AND. (jarray(ij).GE.0) &
+ evolve_mom(ip,ij) = ((parray(ip).GE.0) .AND. (jarray(ij).GE.0)) &
.AND. (parray(ip).LE.pmax) .AND. (jarray(ij).LE.jmax)
ENDDO
ENDDO
CASE('max_degree')
DO ip = 1,local_np+ngp
DO ij = 1, local_nj+ngj
- evolve_mom(ip,ij) = (parray(ip).GE.0) .AND. (jarray(ij).GE.0) &
+ evolve_mom(ip,ij) = ((parray(ip).GE.0) .AND. (jarray(ij).GE.0)) &
.AND. (parray(ip)+2*jarray(ij) .LE. dmax)
ENDDO
ENDDO
CASE DEFAULT
ERROR STOP "closure scheme not recognized (avail: truncation,max_degree)"
END SELECT
-
! Set the nonlinear closure scheme (truncation of sum over n in Sapj)
ALLOCATE(nmaxarray(local_nj))
SELECT CASE(nonlinear_closure)
diff --git a/src/cosolver_interface_mod.F90 b/src/cosolver_interface_mod.F90
index 988c3b0100228343da1ebf6664452d2c982b4167..268058e90d12a0813296d6b3213d4397a0d98c78 100644
--- a/src/cosolver_interface_mod.F90
+++ b/src/cosolver_interface_mod.F90
@@ -1,6 +1,6 @@
module cosolver_interface
! contains the Hermite-Laguerre collision operators solved using COSOlver.
-USE prec_const, ONLY: xp
+USE prec_const, ONLY: xp, mpi_xp_c
IMPLICIT NONE
PRIVATE
PUBLIC :: load_COSOlver_mat, compute_cosolver_coll
@@ -47,11 +47,11 @@ CONTAINS
ENDDO p
IF (num_procs_p .GT. 1) THEN
! Reduce the local_sums to root = 0
- CALL MPI_REDUCE(local_coll, buffer, total_np, MPI_DOUBLE_COMPLEX, MPI_SUM, 0, comm_p, ierr)
+ CALL MPI_REDUCE(local_coll, buffer, total_np, mpi_xp_c, MPI_SUM, 0, comm_p, ierr)
! buffer contains the entire collision term along p, we scatter it between
! the other processes (use of scatterv since Pmax/Np is not an integer)
- CALL MPI_SCATTERV(buffer, rcv_p, dsp_p, MPI_DOUBLE_COMPLEX,&
- TColl_distr, local_np, MPI_DOUBLE_COMPLEX, &
+ CALL MPI_SCATTERV(buffer, rcv_p, dsp_p, mpi_xp_c,&
+ TColl_distr, local_np, mpi_xp_c, &
0, comm_p, ierr)
ELSE
TColl_distr = local_coll
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 594a33c1d105688de72950776fd3177775d01aab..96506661b6cbede2d687034b2878a7b57f8d226a 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -96,15 +96,20 @@ MODULE fourier
REAL(xp), INTENT(IN) :: inv_Nx, inv_Ny
REAL(xp), DIMENSION(local_nky_ptr), INTENT(IN) :: ky_, AA_y
REAL(xp), DIMENSION(local_nkx_ptr), INTENT(IN) :: kx_, AA_x
- COMPLEX(c_xp_c), DIMENSION(local_nky_ptr,local_nkx_ptr), &
- INTENT(IN) :: F_(:,:), G_(:,:)
+ COMPLEX(c_xp_c), DIMENSION(local_nky_ptr,local_nkx_ptr) &
+ :: F_(:,:), G_(:,:)
real(c_xp_r), pointer, INTENT(INOUT) :: sum_real_(:,:)
INTEGER :: ikx,iky
+ !! Anti aliasing
+ DO ikx = 1,local_nkx_ptr
+ F_(:,ikx) = F_(:,ikx)*AA_y(:)*AA_x(ikx)
+ G_(:,ikx) = G_(:,ikx)*AA_y(:)*AA_x(ikx)
+ ENDDO
! First term df/dx x dg/dy
DO ikx = 1,local_nkx_ptr
DO iky = 1,local_nky_ptr
- cmpx_data_f(ikx,iky) = imagu*kx_(ikx)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky)
- cmpx_data_g(ikx,iky) = imagu*ky_(iky)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky)
+ cmpx_data_f(ikx,iky) = imagu*kx_(ikx)*F_(iky,ikx)!*AA_x(ikx)*AA_y(iky)
+ cmpx_data_g(ikx,iky) = imagu*ky_(iky)*G_(iky,ikx)!*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
@@ -120,9 +125,9 @@ MODULE fourier
DO ikx = 1,local_nkx_ptr
DO iky = 1,local_nky_ptr
cmpx_data_f(ikx,iky) = &
- imagu*ky_(iky)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky)
+ imagu*ky_(iky)*F_(iky,ikx)!*AA_x(ikx)*AA_y(iky)
cmpx_data_g(ikx,iky) = &
- imagu*kx_(ikx)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky)
+ imagu*kx_(ikx)*G_(iky,ikx)!*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
#ifdef SINGLE_PRECISION
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 7d6dd5d6cd18779e127938436361972d9c416848..ffcedfdc2081b3eaae89ea28c574be1ad9f7dbd2 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -93,7 +93,7 @@ CONTAINS
END SUBROUTINE geometry_readinputs
subroutine eval_magnetic_geometry
- USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, Ngz,&
+ USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, ngz,&
kxarray, kyarray, set_kparray, nzgrid, zweights_SR, ieven
USE basic, ONLY: speak
USE miller, ONLY: set_miller_parameters, get_miller
@@ -106,7 +106,7 @@ CONTAINS
INTEGER :: eo,iz,iky,ikx
! Allocate arrays
- CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,nzgrid)
+ CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
!
IF( (total_nky .EQ. 1) .AND. (total_nz .EQ. 1)) THEN !1D perp linear run
CALL speak('1D perpendicular geometry')
@@ -142,7 +142,7 @@ CONTAINS
CALL set_kparray(gxx,gxy,gyy,hatB)
DO eo = 1,nzgrid
! Curvature operator (Frei et al. 2022 eq 2.15)
- DO iz = 1,local_nz+Ngz
+ DO iz = 1,local_nz+ngz
G1 = gxx(iz,eo)*gyy(iz,eo)-gxy(iz,eo)*gxy(iz,eo)
G2 = gxx(iz,eo)*gyz(iz,eo)-gxy(iz,eo)*gxz(iz,eo)
G3 = gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo)
@@ -178,7 +178,7 @@ CONTAINS
!
SUBROUTINE eval_salpha_geometry
- USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
+ USE grid, ONLY : local_nz,ngz,zarray,nzgrid
! evaluate s-alpha geometry model
implicit none
REAL(xp) :: z
@@ -186,7 +186,7 @@ CONTAINS
alpha_MHD = 0._xp
DO eo = 1,nzgrid
- DO iz = 1,local_nz+Ngz
+ DO iz = 1,local_nz+ngz
z = zarray(iz,eo)
! metric
@@ -230,14 +230,14 @@ CONTAINS
!
SUBROUTINE eval_zpinch_geometry
- USE grid, ONLY : local_nz,Ngz,zarray,nzgrid
+ USE grid, ONLY : local_nz,ngz,zarray,nzgrid
implicit none
REAL(xp) :: z
INTEGER :: iz, eo
alpha_MHD = 0._xp
DO eo = 1,nzgrid
- DO iz = 1,local_nz+Ngz
+ DO iz = 1,local_nz+ngz
z = zarray(iz,eo)
! metric
@@ -280,13 +280,13 @@ CONTAINS
!--------------------------------------------------------------------------------
! NOT TESTED
subroutine eval_1D_geometry
- USE grid, ONLY : local_nz,Ngz,zarray, nzgrid
+ USE grid, ONLY : local_nz,ngz,zarray, nzgrid
! evaluate 1D perp geometry model
implicit none
REAL(xp) :: z
INTEGER :: iz, eo
DO eo = 1,nzgrid
- DO iz = 1,local_nz+Ngz
+ DO iz = 1,local_nz+ngz
z = zarray(iz,eo)
! metric
@@ -457,31 +457,31 @@ END SUBROUTINE set_ikx_zBC_map
!--------------------------------------------------------------------------------
!
- SUBROUTINE geometry_allocate_mem(local_nky,local_nkx,local_nz,Ngz,nzgrid)
- INTEGER, INTENT(IN) :: local_nky,local_nkx,local_nz,Ngz,nzgrid
+ SUBROUTINE geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
+ INTEGER, INTENT(IN) :: local_nky,local_nkx,local_nz,ngz,nzgrid
! Curvature and geometry
- ALLOCATE( Ckxky(local_nky,local_nkx,local_nz+Ngz,nzgrid))
- ALLOCATE( Jacobian(local_nz+Ngz,nzgrid))
- ALLOCATE( gxx(local_nz+Ngz,nzgrid))
- ALLOCATE( gxy(local_nz+Ngz,nzgrid))
- ALLOCATE( gxz(local_nz+Ngz,nzgrid))
- ALLOCATE( gyy(local_nz+Ngz,nzgrid))
- ALLOCATE( gyz(local_nz+Ngz,nzgrid))
- ALLOCATE( gzz(local_nz+Ngz,nzgrid))
- ALLOCATE( dBdx(local_nz+Ngz,nzgrid))
- ALLOCATE( dBdy(local_nz+Ngz,nzgrid))
- ALLOCATE( dBdz(local_nz+Ngz,nzgrid))
- ALLOCATE( dlnBdz(local_nz+Ngz,nzgrid))
- ALLOCATE( hatB(local_nz+Ngz,nzgrid))
- ! ALLOCATE(Gamma_phipar,(local_nz+Ngz,nzgrid)) (not implemented)
- ALLOCATE( hatR(local_nz+Ngz,nzgrid))
- ALLOCATE( hatZ(local_nz+Ngz,nzgrid))
- ALLOCATE( Rc(local_nz+Ngz,nzgrid))
- ALLOCATE( phic(local_nz+Ngz,nzgrid))
- ALLOCATE( Zc(local_nz+Ngz,nzgrid))
- ALLOCATE( dxdR(local_nz+Ngz,nzgrid))
- ALLOCATE( dxdZ(local_nz+Ngz,nzgrid))
- ALLOCATE(gradz_coeff(local_nz+Ngz,nzgrid))
+ ALLOCATE( Ckxky(local_nky,local_nkx,local_nz+ngz,nzgrid))
+ ALLOCATE( Jacobian(local_nz+ngz,nzgrid))
+ ALLOCATE( gxx(local_nz+ngz,nzgrid))
+ ALLOCATE( gxy(local_nz+ngz,nzgrid))
+ ALLOCATE( gxz(local_nz+ngz,nzgrid))
+ ALLOCATE( gyy(local_nz+ngz,nzgrid))
+ ALLOCATE( gyz(local_nz+ngz,nzgrid))
+ ALLOCATE( gzz(local_nz+ngz,nzgrid))
+ ALLOCATE( dBdx(local_nz+ngz,nzgrid))
+ ALLOCATE( dBdy(local_nz+ngz,nzgrid))
+ ALLOCATE( dBdz(local_nz+ngz,nzgrid))
+ ALLOCATE( dlnBdz(local_nz+ngz,nzgrid))
+ ALLOCATE( hatB(local_nz+ngz,nzgrid))
+ ! ALLOCATE(Gamma_phipar,(local_nz+ngz,nzgrid)) (not implemented)
+ ALLOCATE( hatR(local_nz+ngz,nzgrid))
+ ALLOCATE( hatZ(local_nz+ngz,nzgrid))
+ ALLOCATE( Rc(local_nz+ngz,nzgrid))
+ ALLOCATE( phic(local_nz+ngz,nzgrid))
+ ALLOCATE( Zc(local_nz+ngz,nzgrid))
+ ALLOCATE( dxdR(local_nz+ngz,nzgrid))
+ ALLOCATE( dxdZ(local_nz+ngz,nzgrid))
+ ALLOCATE(gradz_coeff(local_nz+ngz,nzgrid))
END SUBROUTINE geometry_allocate_mem
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index 217b19140e11e77d1896df9c6835910ed693813e..246fd7d3139f213b584b8cb1896f5215d0810827 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -55,26 +55,17 @@ SUBROUTINE update_ghosts_p_mom
USE grid, ONLY: local_na,local_np,local_nj,local_nky,local_nkx,local_nz,&
ngp,ngj,ngz
IMPLICIT NONE
- INTEGER :: ierr, first, last, count
+ INTEGER :: ierr, first, last, count, iz,ikx,iky,ij,igp,ia
+ COMPLEX(xp), DIMENSION(local_na,ngp/2,local_nj+ngj,local_nky,local_nkx,local_nz+ngz) :: buffer_send, buffer_recv
first = 1 + ngp/2
last = local_np + ngp/2
- ! count = local_na*(local_nj+ngj)*local_nky*local_nkx*(local_nz+ngz) ! Number of elements to send
- !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
- ! DO ig = 1,ngp/2
- ! CALL mpi_sendrecv(moments(:,last-(ig-1),:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_R, 14+ig, &
- ! moments(:,first-ig ,:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_L, 14+ig, &
- ! comm0, status, ierr)
- ! ENDDO
- !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
- ! DO ig = 1,ngp/2
- ! CALL mpi_sendrecv(moments(:,first+(ig-1),:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_L, 16+ig, &
- ! moments(:,last + ig ,:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_R, 16+ig, &
- ! comm0, status, ierr)
- ! ENDDO
+
count = (ngp/2)*local_na*(local_nj+ngj)*local_nky*local_nkx*(local_nz+ngz) ! Number of elements to send
+ !!!!!! Send to the right, receive from the left
CALL mpi_sendrecv(moments(:,(last-(ngp/2-1)):(last),:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_R, 14, &
moments(:,(first-ngp/2):(first-1),:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_L, 14, &
comm0, status, ierr)
+ !!!!!!! Send to the left, receive from the right
CALL mpi_sendrecv(moments(:,(first):(first+(ngp/2-1)),:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_L, 16, &
moments(:,(last+1):(last+ngp/2) ,:,:,:,:,updatetlevel), count, mpi_xp_c, nbr_R, 16, &
comm0, status, ierr)
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 495347ebbe17628c42937c6bca34959cc64b1395..f998b3b41f5ce6b203755b283dc0d1a0fe13f8eb 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -6,7 +6,7 @@ MODULE miller
USE basic
USE parallel, ONLY: my_id, num_procs_z, nbr_U, nbr_D, comm0
! use coordinates,only: gcoor, get_dzprimedz
- USE grid, ONLY: local_Nky, local_Nkx, local_nz, Ngz, kyarray, kxarray, zarray, Nz, local_nz_offset
+ USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset
! use discretization
USE lagrange_interpolation
! use par_in, only: beta, sign_Ip_CW, sign_Bt_CW, Npol
@@ -59,12 +59,12 @@ CONTAINS
real(xp), INTENT(INOUT) :: edge_opt ! alpha mhd
real(xp), INTENT(INOUT) :: xpdx_pm_geom ! amplitude mag. eq. pressure grad.
real(xp), INTENT(INOUT) :: C_y, C_xy
- real(xp), dimension(1:local_nz+Ngz,1:2), INTENT(INOUT) :: &
+ real(xp), dimension(local_nz+ngz,nzgrid), INTENT(INOUT) :: &
gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
dBdx_,dBdy_,Bfield_,jacobian_,&
dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
gradz_coeff_
- real(xp), dimension(1:local_Nky,1:local_Nkx,1:local_nz+Ngz,1:2), INTENT(INOUT) :: Ckxky_
+ real(xp), dimension(local_Nky,local_Nkx,local_nz+ngz,nzgrid), INTENT(INOUT) :: Ckxky_
INTEGER :: iz, ikx, iky, eo
! No parameter in gyacomo yet
real(xp) :: sign_Ip_CW=1 ! current sign (only normal current)
@@ -96,8 +96,8 @@ CONTAINS
real(xp), dimension(500*(Npol+1)):: gxx, gxy, gxz, gyy, gyz, gzz, dtheta_dchi_s, dBp_dchi_s, jacobian, dBdx, dBdz
real(xp), dimension(500*(Npol+1)):: g_xx, g_xy, g_xz, g_yy, g_yz, g_zz, tmp_arr_s, dxdR_s, dxdZ_s, K_x, K_y !tmp_arr2
- real(xp), dimension(1:Nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
- real(xp), dimension(1:Nz):: R_out, Z_out, dxdR_out, dxdZ_out
+ real(xp), dimension(1:total_nz):: gxx_out,gxy_out,gxz_out,gyy_out,gyz_out,gzz_out,Bfield_out,jacobian_out, dBdx_out, dBdz_out, chi_out
+ real(xp), dimension(1:total_nz):: R_out, Z_out, dxdR_out, dxdZ_out
real(xp):: d_inv, drPsi, dxPsi, dq_dx, dq_xpsi, R0, Z0, B0, F, D0_full, D1_full, D2_full, D3_full
!real(xp) :: Lnorm, Psi0 ! currently module-wide defined anyway
real(xp):: pprime, ffprime, D0_mid, D1_mid, D2_mid, D3_mid, dx_drho, pi, mu_0, dzprimedz
@@ -435,13 +435,13 @@ CONTAINS
if (edge_opt==0.0) then
!gene z-grid
- chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/Nz,Nz)
+ chi_out=linspace(-pi*Npol,pi*Npol-2*pi*Npol/total_nz,total_nz)
else
!new parallel coordinate chi_out==zprime
!see also tracer_aux.F90
if (Npol>1) ERROR STOP '>> ERROR << Npol>1 has not been implemented for edge_opt=\=0.0'
- do k=1,Nz
- chi_out(k)=sinh((-pi+k*2.*pi/Nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
+ do k=1,total_nz
+ chi_out(k)=sinh((-pi+k*2.*pi/total_nz)*log(edge_opt*pi+sqrt(edge_opt**2*pi**2+1))/pi)/edge_opt
enddo
!transform metrics according to chain rule
do k=1,np_s
@@ -461,73 +461,73 @@ CONTAINS
endif !edge_opt
!interpolate down to GENE z-grid
- call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,Nz)
- call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,Nz)
- call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,Nz)
- call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,Nz)
- call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,Nz)
- call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,Nz)
- call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,Nz)
- call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,Nz)
- call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,Nz)
- call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,Nz)
- call lag3interp(R_s,chi_s,np_s,R_out,chi_out,Nz)
- call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,Nz)
- call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,Nz)
- call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,Nz)
+ call lag3interp(gxx,chi_s,np_s,gxx_out,chi_out,total_nz)
+ call lag3interp(gxy,chi_s,np_s,gxy_out,chi_out,total_nz)
+ call lag3interp(gxz,chi_s,np_s,gxz_out,chi_out,total_nz)
+ call lag3interp(gyy,chi_s,np_s,gyy_out,chi_out,total_nz)
+ call lag3interp(gyz,chi_s,np_s,gyz_out,chi_out,total_nz)
+ call lag3interp(gzz,chi_s,np_s,gzz_out,chi_out,total_nz)
+ call lag3interp(B_s,chi_s,np_s,Bfield_out,chi_out,total_nz)
+ call lag3interp(jacobian,chi_s,np_s,jacobian_out,chi_out,total_nz)
+ call lag3interp(dBdx,chi_s,np_s,dBdx_out,chi_out,total_nz)
+ call lag3interp(dBdz,chi_s,np_s,dBdz_out,chi_out,total_nz)
+ call lag3interp(R_s,chi_s,np_s,R_out,chi_out,total_nz)
+ call lag3interp(Z_s,chi_s,np_s,Z_out,chi_out,total_nz)
+ call lag3interp(dxdR_s,chi_s,np_s,dxdR_out,chi_out,total_nz)
+ call lag3interp(dxdZ_s,chi_s,np_s,dxdZ_out,chi_out,total_nz)
! Fill the interior of the geom arrays with the results
- do eo=1,2
+ do eo=1,nzgrid
DO iz = 1,local_nz
- gxx_(iz+Ngz/2,eo) = gxx_out(iz-local_nz_offset)
- gyy_(iz+Ngz/2,eo) = gyy_out(iz-local_nz_offset)
- gxz_(iz+Ngz/2,eo) = gxz_out(iz-local_nz_offset)
- gyz_(iz+Ngz/2,eo) = gyz_out(iz-local_nz_offset)
- dBdx_(iz+Ngz/2,eo) = dBdx_out(iz-local_nz_offset)
- dBdy_(iz+Ngz/2,eo) = 0.
- gxy_(iz+Ngz/2,eo) = gxy_out(iz-local_nz_offset)
- gzz_(iz+Ngz/2,eo) = gzz_out(iz-local_nz_offset)
- Bfield_(iz+Ngz/2,eo) = Bfield_out(iz-local_nz_offset)
- jacobian_(iz+Ngz/2,eo) = jacobian_out(iz-local_nz_offset)
- dBdz_(iz+Ngz/2,eo) = dBdz_out(iz-local_nz_offset)
- R_hat_(iz+Ngz/2,eo) = R_out(iz-local_nz_offset)
- Z_hat_(iz+Ngz/2,eo) = Z_out(iz-local_nz_offset)
- dxdR_(iz+Ngz/2,eo) = dxdR_out(iz-local_nz_offset)
- dxdZ_(iz+Ngz/2,eo) = dxdZ_out(iz-local_nz_offset)
+ gxx_(iz+ngz/2,eo) = gxx_out(iz+local_nz_offset)
+ gyy_(iz+ngz/2,eo) = gyy_out(iz+local_nz_offset)
+ gxz_(iz+ngz/2,eo) = gxz_out(iz+local_nz_offset)
+ gyz_(iz+ngz/2,eo) = gyz_out(iz+local_nz_offset)
+ dBdx_(iz+ngz/2,eo) = dBdx_out(iz+local_nz_offset)
+ dBdy_(iz+ngz/2,eo) = 0.
+ gxy_(iz+ngz/2,eo) = gxy_out(iz+local_nz_offset)
+ gzz_(iz+ngz/2,eo) = gzz_out(iz+local_nz_offset)
+ Bfield_(iz+ngz/2,eo) = Bfield_out(iz+local_nz_offset)
+ jacobian_(iz+ngz/2,eo) = jacobian_out(iz+local_nz_offset)
+ dBdz_(iz+ngz/2,eo) = dBdz_out(iz+local_nz_offset)
+ R_hat_(iz+ngz/2,eo) = R_out(iz+local_nz_offset)
+ Z_hat_(iz+ngz/2,eo) = Z_out(iz+local_nz_offset)
+ dxdR_(iz+ngz/2,eo) = dxdR_out(iz+local_nz_offset)
+ dxdZ_(iz+ngz/2,eo) = dxdZ_out(iz+local_nz_offset)
ENDDO
- !! UPDATE GHOSTS VALUES (since the miller function in GENE does not)
- CALL update_ghosts_z(gxx_(:,eo))
- CALL update_ghosts_z(gyy_(:,eo))
- CALL update_ghosts_z(gxz_(:,eo))
- CALL update_ghosts_z(gxy_(:,eo))
- CALL update_ghosts_z(gzz_(:,eo))
- CALL update_ghosts_z(Bfield_(:,eo))
- CALL update_ghosts_z(dBdx_(:,eo))
- CALL update_ghosts_z(dBdy_(:,eo))
- CALL update_ghosts_z(dBdz_(:,eo))
- CALL update_ghosts_z(jacobian_(:,eo))
- CALL update_ghosts_z(R_hat_(:,eo))
- CALL update_ghosts_z(Z_hat_(:,eo))
- CALL update_ghosts_z(dxdR_(:,eo))
- CALL update_ghosts_z(dxdZ_(:,eo))
-
- ! Curvature operator (Frei et al. 2022 eq 2.15)
- DO iz = 1,local_nz+Ngz
- G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
- G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
- G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
- Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo)
- Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo)
-
- DO iky = 1,local_Nky
- ky = kyarray(iky)
- DO ikx= 1,local_Nkx
- kx = kxarray(ikx)
- Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)
- ENDDO
+ !! UPDATE GHOSTS VALUES (since the miller function in GENE does not)
+ CALL update_ghosts_z(gxx_(:,eo))
+ CALL update_ghosts_z(gyy_(:,eo))
+ CALL update_ghosts_z(gxz_(:,eo))
+ CALL update_ghosts_z(gxy_(:,eo))
+ CALL update_ghosts_z(gzz_(:,eo))
+ CALL update_ghosts_z(Bfield_(:,eo))
+ CALL update_ghosts_z(dBdx_(:,eo))
+ CALL update_ghosts_z(dBdy_(:,eo))
+ CALL update_ghosts_z(dBdz_(:,eo))
+ CALL update_ghosts_z(jacobian_(:,eo))
+ CALL update_ghosts_z(R_hat_(:,eo))
+ CALL update_ghosts_z(Z_hat_(:,eo))
+ CALL update_ghosts_z(dxdR_(:,eo))
+ CALL update_ghosts_z(dxdZ_(:,eo))
+
+ ! Curvature operator (Frei et al. 2022 eq 2.15)
+ DO iz = 1,local_nz+ngz
+ G1 = gxy_(iz,eo)*gxy_(iz,eo)-gxx_(iz,eo)*gyy_(iz,eo)
+ G2 = gxy_(iz,eo)*gxz_(iz,eo)-gxx_(iz,eo)*gyz_(iz,eo)
+ G3 = gyy_(iz,eo)*gxz_(iz,eo)-gxy_(iz,eo)*gyz_(iz,eo)
+ Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo)
+ Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo)
+
+ DO iky = 1,local_Nky
+ ky = kyarray(iky)
+ DO ikx= 1,local_Nkx
+ kx = kxarray(ikx)
+ Ckxky_(iky, ikx, iz,eo) = Cx*kx + Cy*ky
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff_(iz,eo) = 1._xp / jacobian_(iz,eo) / Bfield_(iz,eo)
ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff_(iz,eo) = 1._xp / jacobian_(iz,eo) / Bfield_(iz,eo)
- ENDDO
ENDDO
contains
@@ -536,12 +536,12 @@ CONTAINS
SUBROUTINE update_ghosts_z(fz_)
IMPLICIT NONE
! INTEGER, INTENT(IN) :: nztot_
- REAL(xp), DIMENSION(1:local_nz+Ngz), INTENT(INOUT) :: fz_
+ REAL(xp), DIMENSION(1:local_nz+ngz), INTENT(INOUT) :: fz_
REAL(xp), DIMENSION(-2:2) :: buff
INTEGER :: status(MPI_STATUS_SIZE), count, last, first
- last = local_nz+Ngz/2
- first= 1 + Ngz/2
- IF(Nz .GT. 1) THEN
+ last = local_nz+ngz/2
+ first= 1 + ngz/2
+ IF(total_nz .GT. 1) THEN
IF (num_procs_z .GT. 1) THEN
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
count = 1 ! one point to exchange
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 7be127dd9e7a782165b4bf824f08bc04465ecaf9..0e306701c7a6000a5514addfe45d1a53aded7d12 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -48,6 +48,7 @@ CONTAINS
ENDIF
IF(Na .EQ. 1) THEN
+ IF(.NOT. ADIAB_E) ERROR STOP "With one species, ADIAB_E must be set to .true. STOP"
IF(my_id.EQ.0) print*, 'Adiabatic electron model -> beta = 0'
beta = 0._xp
ENDIF
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index e2a25dd4696741cc717436f2a87bcada6f186a9c..146e28ab1668e8a44c563d4f9cb915087185cecc 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -3,235 +3,159 @@ MODULE moments_eq_rhs
PUBLIC :: compute_moments_eq_rhs
CONTAINS
-SUBROUTINE compute_moments_eq_rhs
- USE model
- USE array
- USE fields
- 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
- USE basic
- USE closure, ONLY: evolve_mom
- USE prec_const
- USE collision
- USE time_integration
- ! USE species, ONLY: xpdx
- USE geometry, ONLY: gradz_coeff, dlnBdz, Ckxky!, Gamma_phipar
- USE calculus, ONLY: interp_z, grad_z, grad_z2
- ! USE species, ONLY: xpdx
- IMPLICIT NONE
- INTEGER :: ia, iz, iky, ikx, ip ,ij, eo ! counters
- INTEGER :: izi,ipi,iji ! interior points counters
- INTEGER :: p_int, j_int ! loops indices and polynom. degrees
- REAL(xp) :: kx, ky, kperp2
- COMPLEX(xp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
- COMPLEX(xp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
- COMPLEX(xp) :: Ldamp, Fmir
- COMPLEX(xp) :: Mperp, Mpara, Dphi, Dpsi
- COMPLEX(xp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
- COMPLEX(xp) :: i_kx,i_ky
- COMPLEX(xp) :: Napj, RHS, phikykxz, psikykxz
- ! Spatial loops
- 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
- ky = kyarray(iky) ! binormal wavevector
- i_ky = imagu * ky ! binormal derivative
- phikykxz = phi(iky,ikx,izi)
- psikykxz = psi(iky,ikx,izi)
- ! Kinetic loops
- j:DO ij = 1,local_nj ! This loop is from 1 to jmaxi+1
- iji = ij+ngj/2
- j_int = jarray(iji)
- p:DO ip = 1,local_np ! Hermite loop
- 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
- ! Species loop
- a:DO ia = 1,local_na
- Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
- RHS = 0._xp
- IF(evolve_mom(ipi,iji)) THEN ! for the closure scheme
- !! Compute moments_ mixing terms
- ! Perpendicular dynamic
- ! term propto n^{p,j}
- Tnapj = xnapj(ia,ip,ij)* nadiab_moments(ia,ipi, iji, iky,ikx,izi)
- ! term propto n^{p+2,j}
- Tnapp2j = xnapp2j(ia,ip) * nadiab_moments(ia,ipi+pp2,iji, iky,ikx,izi)
- ! term propto n^{p-2,j}
- Tnapm2j = xnapm2j(ia,ip) * nadiab_moments(ia,ipi-pp2,iji, iky,ikx,izi)
- ! term propto n^{p,j+1}
- Tnapjp1 = xnapjp1(ia,ij) * nadiab_moments(ia,ipi, iji+1,iky,ikx,izi)
- ! term propto n^{p,j-1}
- Tnapjm1 = xnapjm1(ia,ij) * nadiab_moments(ia,ipi, iji-1,iky,ikx,izi)
- ! Perpendicular magnetic term (curvature and gradient drifts)
- Mperp = imagu*Ckxky(iky,ikx,izi,eo)&
- *(Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)
- ! Parallel dynamic
- ! ddz derivative for Landau damping term
- Ldamp = xnapp1j(ia,ip) * ddz_napj(ia,ipi+1,iji,iky,ikx,iz) &
- + xnapm1j(ia,ip) * ddz_napj(ia,ipi-1,iji,iky,ikx,iz)
- ! Mirror terms
- Tnapp1j = ynapp1j (ia,ip,ij) * interp_napj(ia,ipi+1,iji ,iky,ikx,iz)
- Tnapp1jm1 = ynapp1jm1(ia,ip,ij) * interp_napj(ia,ipi+1,iji-1,iky,ikx,iz)
- Tnapm1j = ynapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
- Tnapm1jm1 = ynapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
- ! Trapping terms
- Unapm1j = znapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
- Unapm1jp1 = znapm1jp1(ia,ip,ij) * interp_napj(ia,ipi-1,iji+1,iky,ikx,iz)
- Unapm1jm1 = znapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
- ! sum the parallel forces
- Fmir = dlnBdz(izi,eo)*(Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 +&
- Unapm1j + Unapm1jp1 + Unapm1jm1)
- ! Parallel magnetic term (Landau damping and the mirror force)
- Mpara = gradz_coeff(izi,eo)*(Ldamp + Fmir)
- !! Electrical potential term
- IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
- Dphi =i_ky*( xphij(ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
- +xphijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
- +xphijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo) )*phikykxz
- ELSE
- Dphi = 0._xp
- ENDIF
- !! Vector potential term
- IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
- Dpsi =-i_ky*( xpsij (ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
- +xpsijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
- +xpsijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo))*psikykxz
+ SUBROUTINE compute_moments_eq_rhs
+ USE model
+ USE array
+ USE fields
+ 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
+ USE basic
+ USE closure, ONLY: evolve_mom
+ USE prec_const
+ USE collision
+ USE time_integration
+ ! USE species, ONLY: xpdx
+ USE geometry, ONLY: gradz_coeff, dlnBdz, Ckxky!, Gamma_phipar
+ USE calculus, ONLY: interp_z, grad_z, grad_z2
+ ! USE species, ONLY: xpdx
+ IMPLICIT NONE
+ INTEGER :: ia, iz, iky, ikx, ip ,ij, eo ! counters
+ INTEGER :: izi,ipi,iji ! interior points counters
+ INTEGER :: p_int, j_int ! loops indices and polynom. degrees
+ REAL(xp) :: kx, ky, kperp2
+ COMPLEX(xp) :: Tnapj, Tnapp2j, Tnapm2j, Tnapjp1, Tnapjm1 ! Terms from b x gradB and drives
+ COMPLEX(xp) :: Tnapp1j, Tnapm1j, Tnapp1jm1, Tnapm1jm1 ! Terms from mirror force with non adiab moments_
+ COMPLEX(xp) :: Ldamp, Fmir
+ COMPLEX(xp) :: Mperp, Mpara, Dphi, Dpsi
+ COMPLEX(xp) :: Unapm1j, Unapm1jp1, Unapm1jm1 ! Terms from mirror force with adiab moments_
+ COMPLEX(xp) :: i_kx,i_ky
+ COMPLEX(xp) :: Napj, RHS, phikykxz, psikykxz
+ ! Spatial loops
+ 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
+ ky = kyarray(iky) ! binormal wavevector
+ i_ky = imagu * ky ! binormal derivative
+ phikykxz = phi(iky,ikx,izi)
+ psikykxz = psi(iky,ikx,izi)
+ ! Kinetic loops
+ j:DO ij = 1,local_nj ! This loop is from 1 to jmaxi+1
+ iji = ij+ngj/2
+ j_int = jarray(iji)
+ p:DO ip = 1,local_np ! Hermite loop
+ 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
+ ! Species loop
+ a:DO ia = 1,local_na
+ Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
+ RHS = 0._xp
+ IF(evolve_mom(ipi,iji)) THEN ! for the closure scheme
+ !! Compute moments_ mixing terms
+ ! Perpendicular dynamic
+ ! term propto n^{p,j}
+ Tnapj = xnapj(ia,ip,ij)* nadiab_moments(ia,ipi, iji, iky,ikx,izi)
+ ! term propto n^{p+2,j}
+ Tnapp2j = xnapp2j(ia,ip) * nadiab_moments(ia,ipi+pp2,iji, iky,ikx,izi)
+ ! term propto n^{p-2,j}
+ Tnapm2j = xnapm2j(ia,ip) * nadiab_moments(ia,ipi-pp2,iji, iky,ikx,izi)
+ ! term propto n^{p,j+1}
+ Tnapjp1 = xnapjp1(ia,ij) * nadiab_moments(ia,ipi, iji+1,iky,ikx,izi)
+ ! term propto n^{p,j-1}
+ Tnapjm1 = xnapjm1(ia,ij) * nadiab_moments(ia,ipi, iji-1,iky,ikx,izi)
+ ! Perpendicular magnetic term (curvature and gradient drifts)
+ Mperp = imagu*Ckxky(iky,ikx,izi,eo)&
+ *(Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)
+ ! Parallel dynamic
+ ! ddz derivative for Landau damping term
+ Ldamp = xnapp1j(ia,ip) * ddz_napj(ia,ipi+1,iji,iky,ikx,iz) &
+ + xnapm1j(ia,ip) * ddz_napj(ia,ipi-1,iji,iky,ikx,iz)
+ ! Mirror terms
+ Tnapp1j = ynapp1j (ia,ip,ij) * interp_napj(ia,ipi+1,iji ,iky,ikx,iz)
+ Tnapp1jm1 = ynapp1jm1(ia,ip,ij) * interp_napj(ia,ipi+1,iji-1,iky,ikx,iz)
+ Tnapm1j = ynapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
+ Tnapm1jm1 = ynapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
+ ! Trapping terms
+ Unapm1j = znapm1j (ia,ip,ij) * interp_napj(ia,ipi-1,iji ,iky,ikx,iz)
+ Unapm1jp1 = znapm1jp1(ia,ip,ij) * interp_napj(ia,ipi-1,iji+1,iky,ikx,iz)
+ Unapm1jm1 = znapm1jm1(ia,ip,ij) * interp_napj(ia,ipi-1,iji-1,iky,ikx,iz)
+ ! sum the parallel forces
+ Fmir = dlnBdz(izi,eo)*(Tnapp1j + Tnapp1jm1 + Tnapm1j + Tnapm1jm1 +&
+ Unapm1j + Unapm1jp1 + Unapm1jm1)
+ ! Parallel magnetic term (Landau damping and the mirror force)
+ Mpara = gradz_coeff(izi,eo)*(Ldamp + Fmir)
+ !! Electrical potential term
+ IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
+ Dphi =i_ky*( xphij(ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
+ +xphijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+ +xphijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo) )*phikykxz
+ ELSE
+ Dphi = 0._xp
+ ENDIF
+ !! Vector potential term
+ IF ( (p_int .LE. 3) .AND. (p_int .GE. 1) ) THEN ! Kronecker p1 or p3
+ Dpsi =-i_ky*( xpsij (ia,ip,ij)*kernel(ia,iji ,iky,ikx,izi,eo) &
+ +xpsijp1(ia,ip,ij)*kernel(ia,iji+1,iky,ikx,izi,eo) &
+ +xpsijm1(ia,ip,ij)*kernel(ia,iji-1,iky,ikx,izi,eo))*psikykxz
+ ELSE
+ Dpsi = 0._xp
+ ENDIF
+ !! Sum of all RHS terms
+ RHS = &
+ ! Nonlinear term Sapj_ = {phi,f}
+ - Sapj(ia,ip,ij,iky,ikx,iz) &
+ ! Perpendicular magnetic term
+ - Mperp &
+ ! Parallel magnetic term
+ - Mpara &
+ ! Drives (density + temperature gradients)
+ - (Dphi + Dpsi) &
+ ! Collision term
+ + Capj(ia,ip,ij,iky,ikx,iz) &
+ ! Perpendicular pressure effects (electromagnetic term) (TO CHECK)
+ ! - i_ky*beta*xpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
+ ! Parallel drive term (should be negligible, to test)
+ ! !! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
+ ! Numerical perpendicular hyperdiffusion
+ -mu_x*diff_kx_coeff*kx**N_HD*Napj &
+ -mu_y*diff_ky_coeff*ky**N_HD*Napj &
+ ! ! ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
+ -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,iz)
+ !! Velocity space dissipation (should be implemented somewhere else)
+ SELECT CASE(HYP_V)
+ CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
+ IF (p_int .GT. 2) &
+ RHS = RHS - mu_p*diff_p_coeff*p_int**6*Napj
+ IF (j_int .GT. 1) &
+ RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
+ CASE('dvpar4')
+ ! fourth order numerical diffusion in vpar
+ IF(p_int .GE. 4) &
+ ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! (not used often so not parallelized)
+ RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
+ ! + dummy Laguerre diff
+ IF (j_int .GE. 2) &
+ RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
+ CASE DEFAULT
+ END SELECT
ELSE
- Dpsi = 0._xp
+ RHS = 0._xp
ENDIF
- !! Sum of all RHS terms
- RHS = &
- ! Nonlinear term Sapj_ = {phi,f}
- - Sapj(ia,ip,ij,iky,ikx,iz) &
- ! Perpendicular magnetic term
- - Mperp &
- ! Parallel magnetic term
- - Mpara &
- ! Drives (density + temperature gradients)
- - (Dphi + Dpsi) &
- ! Collision term
- + Capj(ia,ip,ij,iky,ikx,iz) &
- ! Perpendicular pressure effects (electromagnetic term) (TO CHECK)
- ! - i_ky*beta*xpdx(ia) * (Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)&
- ! Parallel drive term (should be negligible, to test)
- ! !! -Gamma_phipar(iz,eo)*Tphi*ddz_phi(iky,ikx,iz) &
- ! Numerical perpendicular hyperdiffusion
- -mu_x*diff_kx_coeff*kx**N_HD*Napj &
- -mu_y*diff_ky_coeff*ky**N_HD*Napj &
- ! ! ! Numerical parallel hyperdiffusion "mu_z*ddz**4" see Pueschel 2010 (eq 25)
- -mu_z*diff_dz_coeff*ddzND_napj(ia,ipi,iji,iky,ikx,iz)
- !! Velocity space dissipation (should be implemented somewhere else)
- SELECT CASE(HYP_V)
- CASE('hypcoll') ! GX like Hermite hypercollisions see Mandell et al. 2023 (eq 3.23), unadvised to use it
- IF (p_int .GT. 2) &
- RHS = RHS - mu_p*diff_p_coeff*p_int**6*Napj
- IF (j_int .GT. 1) &
- RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
- CASE('dvpar4')
- ! fourth order numerical diffusion in vpar
- IF(p_int .GE. 4) &
- ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! (not used often so not parallelized)
- RHS = RHS + mu_p*dv4_Hp_coeff(p_int)*moments(ia,ipi-4,iji,iky,ikx,izi,updatetlevel)
- ! + dummy Laguerre diff
- IF (j_int .GE. 2) &
- RHS = RHS - mu_j*diff_j_coeff*j_int**6*Napj
- CASE DEFAULT
- END SELECT
- ELSE
- RHS = 0._xp
- ENDIF
- !! Put RHS in the array
- moments_rhs(ia,ip,ij,iky,ikx,iz,updatetlevel) = RHS
- END DO a
- END DO p
- END DO j
- END DO y
- END DO x
- END DO z
-
- ! print*,'sumabs moments i', SUM(ABS(moments(1,:,:,:,:,:,updatetlevel)))
- ! print*,'moment rhs i 221', moments_rhs(1,1,1,2,2,1,updatetlevel)
- ! ! print*,'sum real nadiabe', SUM(REAL(nadiab_moments(2,:,:,:,:,:)))
- ! print*,'sumreal momrhs i', SUM(REAL(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumimag momrhs i', SUM(IMAG(moments_rhs(1,:,:,:,:,:,:)))
- ! print*,'sumreal phi ', SUM(REAL(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'sumimag phi ', SUM(IMAG(phi(:,:,(1+ngz/2):(local_nz+ngz/2))))
- ! print*,'phi(2,2,1) ', phi(2,2,1+ngz/2)
- ! print*,'sumreal ddznipj ', SUM(REAL(ddz_napj(1,:,:,:,:,:)))
- ! print*,'sumimag ddznipj ', SUM(IMAG(ddz_napj(1,:,:,:,:,:)))
- ! print*,' ddznipj ',(ddz_napj(1,2+ngp/2,2+ngj/2,2,2,1))
- ! print*,' ddzNDnipj ',SUM(REAL(ddzND_Napj(1,:,:,:,:,:)))
- ! print*,'sumreal Capj ', SUM(REAL(Capj(1,:,:,:,:,:)))
- ! print*,'sum phi coeff', SUM(xphij(1,:,:)) + SUM(xphijp1(1,:,:)) + SUM(xphijm1(1,:,:))
- ! print*,'---'
- ! IF(updatetlevel .EQ. 4) stop
- ! stop
-
-END SUBROUTINE compute_moments_eq_rhs
-
-! SUBROUTINE add_Maxwellian_background_terms
-! ! This routine is meant to add the terms rising from the magnetic operator,
-! ! i.e. (B x k_gB) Grad, applied on the background Maxwellian distribution
-! ! (x_a + spar^2)(b x k_gB) GradFaM
-! ! It gives birth to kx=ky=0 sources terms (averages) that hit moments_ 00, 20,
-! ! 40, 01,02, 21 with background gradient dependences.
-! USE prec_const
-! USE time_integration, ONLY : updatetlevel
-! USE species, ONLY: tau_q, k_N, k_T
-! USE array, ONLY: moments_rhs
-! USE grid, ONLY: contains_kx0, contains_ky0, ikx0, iky0,&
-! ia,ias,iae,ip,ips,ipe, ij,ijs,ije, zarray,izs,ize
-! IMPLICIT NONE
-! real(xp), DIMENSION(izs:ize) :: sinz
-!
-! sinz(izs:ize) = SIN(zarray(izs:ize,0))
-!
-! IF(contains_kx0 .AND. contains_ky0) THEN
-! DO ia = ias, iae
-! DO ip = ips,ipe
-! DO ij = ijs,ije
-! SELECT CASE(ij-1)
-! CASE(0) ! j = 0
-! SELECT CASE (ip-1)
-! CASE(0) ! Na00 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * (1.5_xp*k_N(ia) - 1.125_xp*k_T(ia))
-! CASE(2) ! Na20 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * (SQRT2*0.5_xp*k_N(ia) - 2.75_xp*k_T(ia))
-! CASE(4) ! Na40 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * SQRT6*0.75_xp*k_T(ia)
-! END SELECT
-! CASE(1) ! j = 1
-! SELECT CASE (ip-1)
-! CASE(0) ! Na01 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! -tau_q(ia) * sinz(izs:ize) * (k_N(ia) + 3.5_xp*k_T(ia))
-! CASE(2) ! Na21 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! -tau_q(ia) * sinz(izs:ize) * SQRT2*k_T(ia)
-! END SELECT
-! CASE(2) ! j = 2
-! SELECT CASE (ip-1)
-! CASE(0) ! Na02 term
-! moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel) = moments_rhs(ia,ip,ij,iky0,ikx0,izs:ize,updatetlevel)&
-! +tau_q(ia) * sinz(izs:ize) * 2._xp*k_T(ia)
-! END SELECT
-! END SELECT
-! ENDDO
-! ENDDO
-! ENDDO
-! ENDIF
-!
-! END SUBROUTINE
+ !! Put RHS in the array
+ moments_rhs(ia,ip,ij,iky,ikx,iz,updatetlevel) = RHS
+ END DO a
+ END DO p
+ END DO j
+ END DO y
+ END DO x
+ END DO z
+ END SUBROUTINE compute_moments_eq_rhs
END MODULE moments_eq_rhs
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 911f2871a231d8bc518626e581c693181963258d..53b63c26ebb090b4892ae05a1595c22fb9d998bb 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -2,12 +2,13 @@ MODULE processing
USE prec_const, ONLY: xp, imagu, SQRT2, SQRT3, onetwelfth, twothird
USE grid, ONLY: &
local_na, local_np, local_nj, local_nky, local_nkx, local_nz, Ngz,Ngj,Ngp,nzgrid, &
- parray,pmax,ip0, iodd, ieven,&
+ parray,pmax,ip0, iodd, ieven, total_nz, &
CONTAINSp0,ip1,CONTAINSp1,ip2,CONTAINSp2,ip3,CONTAINSp3,&
jarray,jmax,ij0, dmax,&
kyarray, AA_y,&
kxarray, AA_x,&
- zarray, zweights_SR, ieven, iodd, inv_deltaz
+ zarray, zweights_SR, ieven, iodd, inv_deltaz,&
+ kroneck_p0, kroneck_p1
USE fields, ONLY: moments, phi, psi
USE array, ONLY : kernel, nadiab_moments, &
ddz_napj, ddzND_Napj, interp_napj,&
@@ -40,9 +41,10 @@ CONTAINS
ALLOCATE( gflux_x(local_na))
ALLOCATE( hflux_x(local_na))
END SUBROUTINE init_process
+
!------------------------------ HIGH FREQUENCY ROUTINES -------------
! The following routines (nadiab computing, interp and grads) must be
-! as fast as possible since they are called every RK substep.
+! as fast as possible since they are called every RK step
! evaluate the non-adiabatique ion moments
!
! n_{pi} = N^{pj} + kernel(j) /tau_i phi delta_p0
@@ -51,70 +53,34 @@ CONTAINS
IMPLICIT NONE
INTEGER :: ia,ip,ij,iky,ikx,iz
!non adiab moments
- DO iz=1,local_nz+ngz
- DO ikx=1,local_nkx
- DO iky=1,local_nky
- DO ij=1,local_nj+ngj
- DO ip=1,local_np+ngp
- DO ia = 1,local_na
- IF(parray(ip) .EQ. 0) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- + q_tau(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*phi(iky,ikx,iz)
- ELSEIF( (parray(ip) .EQ. 1) .AND. (beta .GT. 0) ) THEN
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel) &
- - q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,iky,ikx,iz,ieven)*psi(iky,ikx,iz)
- ELSE
- nadiab_moments(ia,ip,ij,iky,ikx,iz) = moments(ia,ip,ij,iky,ikx,iz,updatetlevel)
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
- ENDDO
+ ! default : same as moments
+ nadiab_moments(:,:,:,:,:,:) = moments(:,:,:,:,:,:,updatetlevel)
+ ! add phi and psi contribution
+ IF (CONTAINSp0) THEN
+ DO ij=1,local_nj+ngj
+ DO ia = 1,local_na
+ nadiab_moments(ia,ip0,ij,:,:,:) = moments(ia,ip0,ij,:,:,:,updatetlevel) &
+ +q_tau(ia)*kernel(ia,ij,:,:,:,ieven)*phi(:,:,:)
+ ENDDO
+ ENDDO
+ ENDIF
+ IF (CONTAINSp1 .AND. (beta .GT. 0)) THEN
+ DO ij=1,local_nj+ngj
+ DO ia = 1,local_na
+ nadiab_moments(ia,ip1,ij,:,:,:) = moments(ia,ip1,ij,:,:,:,updatetlevel) &
+ -q_o_sqrt_tau_sigma(ia)*kernel(ia,ij,:,:,:,iodd)*psi(:,:,:)
+ ENDDO
+ ENDDO
+ ENDIF
END SUBROUTINE compute_nadiab_moments
- ! z grid gradients
- ! SUBROUTINE compute_gradients_z
- ! IMPLICIT NONE
- ! INTEGER :: eo, p_int, j_int, ia,ip,ij,iky,ikx,iz,izi
- ! COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
- ! COMPLEX(xp), DIMENSION(local_nz) :: f_out
- ! ! Compute z first derivative
- ! DO iz=1,local_nz+ngz
- ! izi = iz+ngz/2
- ! DO ikx=1,local_nkx
- ! DO iky=1,local_nky
- ! DO ij=1,local_nj+ngj
- ! DO ip=1,local_np+ngp
- ! DO ia = 1,local_na
- ! ddz_napj(ia,ip,ij,iky,ikx,iz) = inv_deltaz *(&
- ! +onetwelfth*nadiab_moments(ia,ip,ij,iky,ikx,izi-2)&
- ! -twothird*nadiab_moments(ia,ip,ij,iky,ikx,izi-1)&
- ! +twothird*nadiab_moments(ia,ip,ij,iky,ikx,izi+1)&
- ! -onetwelfth*nadiab_moments(ia,ip,ij,iky,ikx,izi-2)&
- ! )
- ! ddzND_Napj(ia,ip,ij,iky,ikx,iz) = inv_deltaz**4 *(&
- ! +1._xp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
- ! -4._xp*moments(ia,ip,ij,iky,ikx,izi-1,updatetlevel)&
- ! +6._xp*moments(ia,ip,ij,iky,ikx,izi ,updatetlevel)&
- ! -4._xp*moments(ia,ip,ij,iky,ikx,izi+1,updatetlevel)&
- ! +1._xp*moments(ia,ip,ij,iky,ikx,izi-2,updatetlevel)&
- ! )
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! ENDDO
- ! END SUBROUTINE compute_gradients_z
-
! ! z grid gradients
SUBROUTINE compute_gradients_z
IMPLICIT NONE
INTEGER :: eo, p_int, ia,ip,ij,iky,ikx,iz
COMPLEX(xp), DIMENSION(local_nz+ngz) :: f_in
COMPLEX(xp), DIMENSION(local_nz) :: f_out
+ IF(total_nz .GT. 4) THEN
DO ikx = 1,local_nkx
DO iky = 1,local_nky
DO ij = 1,local_nj+ngj
@@ -129,23 +95,22 @@ CONTAINS
! Compute z first derivative
f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
CALL grad_z(eo,local_nz,ngz,inv_deltaz,f_in,f_out)
- ddz_napj(ia,ip,ij,iky,ikx,:) = f_out(:)
+ ddz_napj(ia,ip,ij,iky,ikx,:) = f_out
! Parallel numerical diffusion
- IF (HDz_h) THEN
- f_in = nadiab_moments(ia,ip,ij,iky,ikx,:)
- ELSE
+ IF (.NOT. HDz_h) &
f_in = moments(ia,ip,ij,iky,ikx,:,updatetlevel)
- ENDIF
CALL grad_z4(local_nz,ngz,inv_deltaz,f_in,f_out)
! get output
- DO iz = 1,local_nz
- ddzND_Napj(ia,ip,ij,iky,ikx,iz) = f_out(iz)
- ENDDO
+ ddzND_Napj(ia,ip,ij,iky,ikx,:) = f_out
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
+ ELSE ! 2D Z-pinch
+ ddz_napj = 0._xp
+ ddzND_Napj = 0._xp
+ ENDIF
END SUBROUTINE compute_gradients_z
! z grid interpolation