diff --git a/src/auxval.F90 b/src/auxval.F90
index 5e69b5ebfaaea198ad9a732379c32cdc0af50f40..eec635940a18e18a89b106a238e82096dd840a7f 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -41,7 +41,7 @@ subroutine auxval
CALL build_dv4Hp_table
! set the closure scheme in use
CALL set_closure_model
-
+
#ifdef TEST_SVD
! If we want to test SVD decomposition etc.
CALL init_CLA(local_nky,local_np*local_nj)
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index 299fc626c1bebb700a5fdccbac4415c2e7b3356d..f89bb1f4350a57cb5d18c1ee3dab57cd3838c498 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -54,13 +54,16 @@ SUBROUTINE set_closure_model
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)) &
- .AND. (parray(ip)+2*jarray(ij) .LE. dmax)
+ evolve_mom(ip,ij) = ((parray(ip).GE.0) .AND. (jarray(ij).GE.0)) .AND. (parray(ip)+2*jarray(ij) .LE. dmax)
+ ! evolve_mom(ip,ij) = (parray(ip).GE.0) .AND. (jarray(ij).GE.0)
+ ! evolve_mom(ip,ij) = ((parray(ip)+2*jarray(ij)) .LE. dmax)
+ ! evolve_mom(ip,ij) = evolve_mom(ip,ij) .AND. ((parray(ip)+2*jarray(ij)) .LE. dmax)
ENDDO
ENDDO
CASE DEFAULT
ERROR STOP "closure scheme not recognized (avail: truncation,max_degree,monomial)"
END SELECT
+
! If monomial truncation, setup the coefficients required
SELECT CASE(hierarchy_closure)
CASE('monomial')
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index d5aca96008b8b22df9ee567a315eadf2cccd45eb..13f7a1020062958db17b28f48526257a38383083 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -97,16 +97,19 @@ SUBROUTINE diagnose_full(kstep)
cstep,iframe0d,iframe3d,iframe5d,crashed
USE grid, ONLY: &
parray_full,pmax,jarray_full,jmax,&
- kyarray_full,kxarray_full,zarray_full
+ kyarray_full,kxarray_full,zarray_full, ngz, total_nz, local_nz, ieven
+ USE geometry, ONLY: gxx, gxy, gxz, gyy, gyz, gzz, &
+ hatR, hatZ, hatB, dBdx, dBdy, dBdz, Jacobian, gradz_coeff
USE diagnostics_par
USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf!, putarrnd ! Routine de merde, jamais l'utiliser
USE array
USE model, ONLY: EM
- USE parallel, ONLY: my_id, comm0
+ USE parallel, ONLY: my_id, comm0, gather_z
USE collision, ONLY: coll_outputinputs
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
INTEGER, parameter :: BUFSIZE = 2
+ REAL(xp), DIMENSION(total_nz) :: Az_full ! full z array for metric output
INTEGER :: rank = 0, ierr
INTEGER :: dims(1) = (/0/)
!____________________________________________________________________________
@@ -140,20 +143,34 @@ SUBROUTINE diagnose_full(kstep)
CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
! Metric info
CALL creatg(fidres, "/data/metric", "Metric data")
- ! CALL putarrnd(fidres, "/data/metric/gxx", gxx((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gxy", gxy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gxz", gxz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gyy", gyy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gyz", gyz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gzz", gzz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/hatR", hatR((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/hatZ", hatZ((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/hatB", hatB((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/dBdx", dBdx((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/dBdy", dBdy((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/dBdz", dBdz((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
- ! CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff((1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 1/))
+ CALL gather_z(gxx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxx", Az_full, "gxx", ionode =0)
+ CALL gather_z(gxy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxy", Az_full, "gxy", ionode =0)
+ CALL gather_z(gxz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gxz", Az_full, "gxz", ionode =0)
+ CALL gather_z(gyy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gyy", Az_full, "gyy", ionode =0)
+ CALL gather_z(gyz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gyz", Az_full, "gyz", ionode =0)
+ CALL gather_z(gzz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gzz", Az_full, "gzz", ionode =0)
+ CALL gather_z(hatR((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatR", Az_full, "hatR", ionode =0)
+ CALL gather_z(hatZ((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatZ", Az_full, "hatZ", ionode =0)
+ CALL gather_z(hatB((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/hatB", Az_full, "hatB", ionode =0)
+ CALL gather_z(dBdx((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdx", Az_full, "dBdx", ionode =0)
+ CALL gather_z(dBdy((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdy", Az_full, "dBdy", ionode =0)
+ CALL gather_z(dBdz((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/dBdz", Az_full, "dBdz", ionode =0)
+ CALL gather_z(Jacobian((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/Jacobian", Az_full, "Jacobian", ionode =0)
+ CALL gather_z(gradz_coeff((1+ngz/2):(local_nz+ngz/2),ieven),Az_full,local_nz,total_nz)
+ CALL putarr(fidres, "/data/metric/gradz_coeff", Az_full, "gradz_coeff", ionode =0)
! CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),:), (/1, 1, 3/))
! CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,(1+ngj/2):(local_nj+ngj/2),1:local_nky,1:local_nkx,(1+ngz/2):(local_nz+ngz/2),1), (/1, 2, 4/))
! var0d group (gyro transport)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index f17e38be100c0262dba2d8e6b982e37cbd5d37c8..1c93c6601ed130dfd9451a4a2a358b1cdae97774 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -10,7 +10,7 @@ implicit none
REAL(xp), PUBLIC, PROTECTED :: q0 = 1.4_xp ! safety factor
REAL(xp), PUBLIC, PROTECTED :: shear = 0._xp ! magnetic field shear
REAL(xp), PUBLIC, PROTECTED :: eps = 0.18_xp ! inverse aspect ratio
- REAL(xp), PUBLIC, PROTECTED :: alpha_MHD = 0 ! shafranov shift effect alpha = -q2 R dbeta/dr
+ REAL(xp), PUBLIC, PROTECTED :: alpha_MHD = 0._xp ! shafranov shift effect alpha = -q2 R dbeta/dr
! parameters for Miller geometry
REAL(xp), PUBLIC, PROTECTED :: kappa = 1._xp ! elongation (1 for circular)
REAL(xp), PUBLIC, PROTECTED :: s_kappa = 0._xp ! r normalized derivative skappa = r/kappa dkappa/dr
@@ -21,6 +21,7 @@ implicit none
! to apply shift in the parallel z-BC if shearless
REAL(xp), PUBLIC, PROTECTED :: shift_y = 0._xp ! for Arno <3
REAL(xp), PUBLIC, PROTECTED :: Npol = 1._xp ! number of poloidal turns (real for 3D zpinch studies)
+ LOGICAL , PUBLIC, PROTECTED :: PB_PHASE = .false. ! To activate the parallel boundary phase factor
! Chooses the type of parallel BC we use for the unconnected kx modes (active for non-zero shear only)
! 'periodic' : Connect a disconnected kx to a mode on the other cadran
! 'dirichlet' : Connect a disconnected kx to 0
@@ -31,9 +32,9 @@ implicit none
! GENE unused additional parameters for miller_mod
REAL(xp), PUBLIC, PROTECTED :: edge_opt = 0._xp ! meant to redistribute the points in z
- REAL(xp), PUBLIC, PROTECTED :: major_R = 1._xp ! major radius
+ REAL(xp), PUBLIC, PROTECTED :: major_R = 1.0_xp ! major radius
REAL(xp), PUBLIC, PROTECTED :: major_Z = 0._xp ! vertical elevation
- REAL(xp), PUBLIC, PROTECTED :: xpdx_pm_geom = 0._xp ! amplitude mag. eq. pressure grad.
+ REAL(xp), PUBLIC, PROTECTED :: xpdx_pm_geom = 1._xp ! amplitude mag. eq. pressure grad.
REAL(xp), PUBLIC, PROTECTED :: C_y = 0._xp ! defines y coordinate : Cy (q theta - phi)
REAL(xp), PUBLIC, PROTECTED :: C_xy = 1._xp ! defines x coordinate : B = Cxy Vx x Vy
@@ -44,7 +45,7 @@ implicit none
! Jacobian
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
COMPLEX(xp), PUBLIC, PROTECTED :: iInt_Jacobian ! Inverse integrated Jacobian
- ! Metric
+ ! Metric (local arrays)
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
@@ -55,6 +56,10 @@ implicit none
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
REAL(xp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+ ! full arrays for output
+ REAL(xp), PUBLIC, DIMENSION(:), ALLOCATABLE :: gxx_full, gxy_full, gxz_full, gyy_full, gyz_full, gzz_full
+ REAL(xp), PUBLIC, DIMENSION(:), ALLOCATABLE :: dxdr_full, dxdZ_full, Rc_full, phic_full, Zc_full
+ REAL(xp), PUBLIC, DIMENSION(:), ALLOCATABLE :: dBdx_full, dBdy_full, dBdz_full, dlnBdz_full, hatB_full, hatR_full, hatZ_full, gradz_coeff_full
! Array to map the index of mode (kx,ky,-pi) to (kx+2pi*s*ky,ky,pi) for sheared periodic boundary condition
INTEGER, PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ikx_zBC_L, ikx_zBC_R
! Geometric factor in front of the parallel phi derivative (not implemented)
@@ -76,7 +81,7 @@ CONTAINS
IMPLICIT NONE
NAMELIST /GEOMETRY/ geom, q0, shear, eps,&
kappa, s_kappa,delta, s_delta, zeta, s_zeta,& ! For miller
- parallel_bc, shift_y, Npol
+ parallel_bc, shift_y, Npol, PB_PHASE
READ(lu_in,geometry)
IF(shear .NE. 0._xp) SHEARED = .true.
SELECT CASE(parallel_bc)
@@ -98,6 +103,8 @@ CONTAINS
USE basic, ONLY: speak
USE miller, ONLY: set_miller_parameters, get_miller
USE calculus, ONLY: simpson_rule_z
+ USE model, ONLY: beta
+ USE species, ONLY: Ptot
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(xp) :: kx,ky
@@ -107,6 +114,10 @@ CONTAINS
! Allocate arrays
CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
+
+ ! Set MHD pressure coefficient, flagged by model module's MHD_PD
+ alpha_MHD = q0**2*beta*Ptot
+
!
IF( (total_nky .EQ. 1) .AND. (total_nz .EQ. 1)) THEN !1D perp linear run
CALL speak('1D perpendicular geometry')
@@ -114,8 +125,9 @@ CONTAINS
ELSE
SELECT CASE(geom)
CASE('s-alpha')
- IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for s-alpha geometry"
CALL speak('s-alpha geometry')
+ IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for s-alpha geometry"
+ IF(MODULO(FLOOR(Npol),2) .EQ. 0) ERROR STOP "Npol must be odd for s-alpha"
call eval_salpha_geometry
CASE('Z-pinch','z-pinch','Zpinch','zpinch')
CALL speak('Z-pinch geometry')
@@ -128,6 +140,7 @@ CONTAINS
CASE('miller')
CALL speak('Miller geometry')
IF(FLOOR(Npol) .NE. CEILING(Npol)) ERROR STOP "ERROR STOP: Npol must be integer for Miller geometry"
+ IF(MODULO(FLOOR(Npol),2) .EQ. 0) ERROR STOP "Npol must be odd for Miller"
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
call get_miller(eps,major_R,major_Z,q0,shear,FLOOR(Npol),alpha_MHD,edge_opt,&
C_y,C_xy,xpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
@@ -185,7 +198,6 @@ CONTAINS
implicit none
REAL(xp) :: z
INTEGER :: iz, eo
- alpha_MHD = 0._xp
DO eo = 1,nzgrid
DO iz = 1,local_nz+ngz
@@ -236,7 +248,6 @@ CONTAINS
implicit none
REAL(xp) :: z
INTEGER :: iz, eo
- alpha_MHD = 0._xp
DO eo = 1,nzgrid
DO iz = 1,local_nz+ngz
@@ -315,10 +326,10 @@ CONTAINS
SUBROUTINE set_ikx_zBC_map
USE grid, ONLY: local_nky,total_nkx,contains_zmin,contains_zmax, Nexc,&
- local_nky_offset
+ local_nky_offset, kx_max, kx_min, kyarray, kxarray
USE prec_const, ONLY: imagu, pi
IMPLICIT NONE
- ! REAL(xp) :: shift
+ REAL(xp) :: shift
INTEGER :: ikx,iky, mn_y
ALLOCATE(ikx_zBC_L(local_nky,total_nkx))
ALLOCATE(ikx_zBC_R(local_nky,total_nkx))
@@ -348,13 +359,18 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol turns
- ! shift = 2._xp*PI*shear*kyarray(iky)*Npol
+ shift = 2._xp*PI*shear*kyarray(iky)*Npol
DO ikx = 1,total_nkx
! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc
ikx_zBC_L(iky,ikx) = ikx-mn_y*Nexc
+ ! 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
- ! IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN
- IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN ! outside of the frequ 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)
+ ! IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN
SELECT CASE(parallel_bc)
CASE ('dirichlet')! connected to 0
ikx_zBC_L(iky,ikx) = -99
@@ -367,7 +383,8 @@ CONTAINS
ENDDO
! phase present in GENE from a shift of the x origin by Lx/2 (useless?)
! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_L(iky) = (-1._xp)**(Nexc*mn_y)*EXP(imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
+ IF (PB_PHASE) &
+ pb_phase_L(iky) = (-1._xp)**(Nexc*mn_y-1)*EXP(imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
@@ -378,13 +395,15 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol
- ! shift = 2._xp*PI*shear*kyarray(iky)*Npol
+ shift = 2._xp*PI*shear*kyarray(iky)*Npol
DO ikx = 1,total_nkx
! Usual formula for shifting indices
ikx_zBC_R(iky,ikx) = ikx+mn_y*Nexc
+ ! 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( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain
+ IF( (kxarray(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') ! connected to 0
ikx_zBC_R(iky,ikx) = -99
@@ -398,7 +417,8 @@ CONTAINS
ENDDO
! phase present in GENE from a shift ofthe x origin by Lx/2 (useless?)
! We also put the user defined shift in the y direction (see Volcokas et al. 2022)
- pb_phase_R(iky) = (-1._xp)**(Nexc*mn_y)*EXP(-imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
+ IF (PB_PHASE) &
+ pb_phase_R(iky) = (-1._xp)**(Nexc*mn_y-1)*EXP(-imagu*REAL(mn_y,xp)*2._xp*pi*shift_y)
ENDDO
ENDIF
! Option for disconnecting every modes, viz. connecting all boundary to 0
@@ -407,12 +427,14 @@ CONTAINS
! write(*,*) kxarray
! write(*,*) kyarray
! write(*,*) 'ikx_zBC_L :-----------'
- ! DO iky = ikys,ikye
+ ! DO iky = 1,local_nky
+ ! print*, 'iky=', iky
! print*, ikx_zBC_L(iky,:)
! enddo
! print*, pb_phase_L
! write(*,*) 'ikx_zBC_R :-----------'
- ! DO iky = ikys,ikye
+ ! DO iky = 1,local_nky
+ ! print*, 'iky=', iky
! print*, ikx_zBC_R(iky,:)
! enddo
! print*, pb_phase_R
@@ -475,7 +497,6 @@ END SUBROUTINE set_ikx_zBC_map
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))
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 89f533248a7f2d6624f7c2263d5ccc5fd65360f1..6aa9ad31d44793bff3df6d21586d3e181546db07 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -438,13 +438,13 @@ CONTAINS
ELSE ! Build apprpopriate grid
deltakx = 2._xp*PI/Lx
IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
- kx_max = (total_nkx/2)*deltakx
- kx_min = -kx_max+deltakx
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
DO ikx = 1,total_nkx
kxarray_full(ikx) = deltakx*REAL(MODULO(ikx-1,total_nkx/2)-(total_nkx/2)*FLOOR(2.*real(ikx-1)/real(total_nkx)),xp)
IF (ikx .EQ. total_nkx/2+1) kxarray_full(ikx) = -kxarray_full(ikx)
END DO
+ kx_max = MAXVAL(kxarray_full)!(total_nkx/2)*deltakx
+ 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
@@ -461,7 +461,7 @@ CONTAINS
ENDIF
END DO
ELSE ! Odd number of kx (-2 -1 0 1 2)
- kx_max = (total_nkx-1)/2*deltakx
+ ERROR STOP "Gyacomo is safer with a even Kx number"
ENDIF
ENDIF
! Orszag 2/3 filter
diff --git a/src/inital.F90 b/src/inital.F90
index 3dd50a2a91d73c53b78b5c1e82455f2a47a79aea..ec4719370c28c41e535a79e5a9270035ded9398e 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -466,11 +466,8 @@ SUBROUTINE init_ricci
USE prec_const, ONLY: xp, imagu
USE initial_par,ONLY: iseed, init_noiselvl, init_background
USE model, ONLY: LINEARITY
- USE parallel, ONLY: my_id
IMPLICIT NONE
-
- COMPLEX(xp), DIMENSION(186,52) :: ricci_mat_real, ricci_mat_imag, ricci_face
- REAL(xp) :: tmp_real, tmp_imag
+ COMPLEX(xp), DIMENSION(186,52) :: ricci_mat_real, ricci_mat_imag
COMPLEX(xp) :: scaling
INTEGER :: ia,ip,ij,ikx,iky,iz
! open data file
@@ -505,7 +502,7 @@ SUBROUTINE init_ricci
END DO
END DO
END DO
- print*, sum(moments)
+
IF ( contains_ky0 ) THEN
DO ip=1+ngp/2,local_np+ngp/2
DO ij=1+ngj/2,local_nj+ngj/2
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index f998b3b41f5ce6b203755b283dc0d1a0fe13f8eb..f6b9959fb8882dca4f07e209e9f8a7c9b63b65b5 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -1,4 +1,4 @@
-!! This source has been adapted from GENE https://genecode.org/ !!
+!! This source has been adapted from the GENE code https://genecode.org/ !!
!>Implementation of the local equilibrium model of [R.L. Miller et al., PoP 5, 973 (1998)
!>and [J. Candy, PPCF 51, 105009 (2009)]
MODULE miller
@@ -56,7 +56,7 @@ CONTAINS
real(xp), INTENT(INOUT) :: shat ! safetyfactor
INTEGER, INTENT(IN) :: Npol ! number of poloidal turns
real(xp), INTENT(INOUT) :: amhd ! alpha mhd
- real(xp), INTENT(INOUT) :: edge_opt ! alpha mhd
+ real(xp), INTENT(INOUT) :: edge_opt ! optimization of point placement
real(xp), INTENT(INOUT) :: xpdx_pm_geom ! amplitude mag. eq. pressure grad.
real(xp), INTENT(INOUT) :: C_y, C_xy
real(xp), dimension(local_nz+ngz,nzgrid), INTENT(INOUT) :: &
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 901c7fe4ca3601397a26a8ab1c86b8a9691bfd41..fffedaba7f6315ae60c9e6fdbecb57cc97659b4a 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -23,9 +23,11 @@ MODULE model
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
- REAL(xp), PUBLIC, PROTECTED :: tau_e = 1.0 ! electron temperature ratio for adiabatic electrons
+ REAL(xp), PUBLIC, PROTECTED :: tau_e = 1.0 ! electron temperature ratio for adiabatic electrons
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .false. ! Electromagnetic effects flag
+ LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .false. ! MHD pressure drift
+
! Removes Landau damping in temperature and higher equation (Ivanov 2022)
LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false.
PUBLIC :: model_readinputs, model_outputinputs
@@ -34,14 +36,14 @@ CONTAINS
SUBROUTINE model_readinputs
! Read the input parameters
- USE basic, ONLY: lu_in
- USE parallel, ONLY: my_id,num_procs_p
+ USE basic, ONLY: lu_in, speak
+ USE parallel, ONLY: num_procs_p
USE prec_const
IMPLICIT NONE
NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, &
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, tau_e
+ nu, k_gB, k_cB, lambdaD, MHD_PD, beta, ADIAB_E, tau_e
READ(lu_in,model_par)
@@ -51,18 +53,17 @@ CONTAINS
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'
+ CALL speak('Adiabatic electron model -> beta = 0')
beta = 0._xp
ENDIF
IF(beta .GT. 0) THEN
- IF(my_id.EQ.0) print*, 'Electromagnetic effects are included'
+ CALL speak('Electromagnetic effects are included')
EM = .TRUE.
ENDIF
END SUBROUTINE model_readinputs
-
SUBROUTINE model_outputinputs(fid)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach, creatd
@@ -87,6 +88,7 @@ CONTAINS
CALL attach(fid, TRIM(str), "k_gB", k_gB)
CALL attach(fid, TRIM(str), "k_cB", k_cB)
CALL attach(fid, TRIM(str), "lambdaD", lambdaD)
+ CALL attach(fid, TRIM(str), "MHD_PD", MHD_PD)
CALL attach(fid, TRIM(str), "beta", beta)
CALL attach(fid, TRIM(str), "ADIAB_E", ADIAB_E)
CALL attach(fid, TRIM(str), "tau_e", tau_e)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index f7c837708c2d5e9dc7df0a3ae996632fd3ba1084..8e29e60c3c31cef6b0bfb2ff6b3769cc25a80213 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -16,10 +16,9 @@ CONTAINS
USE prec_const
USE collision
USE time_integration
- ! USE species, ONLY: xpdx
+ USE species, ONLY: Ptot
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
@@ -69,8 +68,8 @@ CONTAINS
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)&
+ ! Perpendicular magnetic term (curvature, gradient drifts and alpha MHD pressure drift)
+ Mperp = imagu*(Ckxky(iky,ikx,izi,eo) - beta*Ptot*ky)&
*(Tnapj + Tnapp2j + Tnapm2j + Tnapjp1 + Tnapjm1)
! Parallel dynamic
! ddz derivative for Landau damping term
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index 0739c3a2a9ab5b8db638e57f9fbc92c5d384aa5e..5f7d306ab7281fb4a75e508f3bec3cf796e6293c 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -1,5 +1,5 @@
MODULE parallel
- use prec_const, ONLY : xp, mpi_xp_c
+ use prec_const, ONLY : xp, mpi_xp_c, mpi_xp_r
USE mpi
IMPLICIT NONE
! Auxiliary variables
@@ -30,6 +30,7 @@ MODULE parallel
! recieve and displacement counts for gatherv
INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_p, dsp_p
INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_y, dsp_y
+ INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_z, dsp_z
INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zy, dsp_zy
INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_zp, dsp_zp
INTEGER, DIMENSION(:), ALLOCATABLE :: rcv_yp, dsp_yp
@@ -120,12 +121,20 @@ CONTAINS
dsp_p(i_) =dsp_p(i_-1) + rcv_p(i_-1)
END DO
!!!!!! XYZ gather variables
+ ! y reduction at constant x,y,z,j
ALLOCATE(rcv_y(num_procs_ky),dsp_y(num_procs_ky)) !Displacement sizes for balance diagnostic
CALL MPI_ALLGATHER(nky_loc,1,MPI_INTEGER,rcv_y,1,MPI_INTEGER,comm_ky,ierr)
dsp_y(1)=0
DO i_=2,num_procs_ky
dsp_y(i_) =dsp_y(i_-1) + rcv_y(i_-1)
END DO
+ ! z reduction at constant x,y,z,j
+ ALLOCATE(rcv_z(num_procs_z),dsp_z(num_procs_z)) !Displacement sizes for balance diagnostic
+ CALL MPI_ALLGATHER(nz_loc,1,MPI_INTEGER,rcv_z,1,MPI_INTEGER,comm_z,ierr)
+ dsp_z(1)=0
+ DO i_=2,num_procs_z
+ dsp_z(i_) =dsp_z(i_-1) + rcv_z(i_-1)
+ END DO
!! Z reduction for full slices of y data but constant x
ALLOCATE(rcv_zy(num_procs_z),dsp_zy(num_procs_z)) !Displacement sizes for balance diagnostic
CALL MPI_ALLGATHER(nz_loc*nky_tot,1,MPI_INTEGER,rcv_zy,1,MPI_INTEGER,comm_z,ierr)
@@ -173,6 +182,31 @@ CONTAINS
CALL attach(fid, TRIM(str), "Np_z", num_procs_z)
END SUBROUTINE parallel_ouptutinputs
+ !!!! Gather a field in z coordinates on rank 0 !!!!!
+ SUBROUTINE gather_z(field_loc,field_tot,nz_loc,nz_tot)
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: nz_loc,nz_tot
+ REAL(xp), DIMENSION(nz_loc), INTENT(IN) :: field_loc
+ REAL(xp), DIMENSION(nz_tot), INTENT(OUT) :: field_tot
+ REAL(xp), DIMENSION(nz_loc) :: buffer_zl !full y, local z
+ REAL(xp), DIMENSION(nz_tot) :: buffer_zt !full z
+ INTEGER :: snd_z, root_z
+
+ snd_z = nz_loc ! Number of points to send along z (full y)
+ root_z = 0
+ buffer_zl = field_loc
+ if(num_procs_z .GT. 1) THEN
+ CALL MPI_GATHERV(buffer_zl, snd_z, mpi_xp_r, &
+ buffer_zt, rcv_z, dsp_z, mpi_xp_r, &
+ root_z, comm_z, ierr)
+ ! ID 0 (the one who output) rebuild the whole array
+ IF(my_id .EQ. 0) &
+ field_tot(1:nz_tot) = buffer_zt(1:nz_tot)
+ ELSE
+ field_tot = field_loc
+ ENDIF
+ END SUBROUTINE gather_z
+
!!!! Gather a field in spatial coordinates on rank 0 !!!!!
SUBROUTINE gather_xyz(field_loc,field_tot,nky_loc,nky_tot,nkx_tot,nz_loc,nz_tot)
IMPLICIT NONE
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index ae3bfdf31103400741e3fbf4e7730562b6e0c384..ce1efa47661cc1a0c5cb03aeaf34bf150b640350 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -24,13 +24,14 @@ MODULE prec_const
INTEGER, PARAMETER :: xp = REAL32
INTEGER, PARAMETER :: c_xp_c = C_FLOAT_COMPLEX
INTEGER, PARAMETER :: c_xp_r = C_FLOAT
+ INTEGER, PARAMETER :: mpi_xp_r = MPI_FLOAT
INTEGER, PARAMETER :: mpi_xp_c = MPI_COMPLEX
#else
INTEGER, PARAMETER :: xp = REAL64
INTEGER, PARAMETER :: c_xp_c = C_DOUBLE_COMPLEX
INTEGER, PARAMETER :: c_xp_r = C_DOUBLE
+ INTEGER, PARAMETER :: mpi_xp_r = MPI_DOUBLE
INTEGER, PARAMETER :: mpi_xp_c = MPI_DOUBLE_COMPLEX
-
#endif
! Auxiliary variables (unused)
INTEGER, private :: xp_r, xp_p !Range and precision of single
diff --git a/src/species_mod.F90 b/src/species_mod.F90
index 071ea9101a591b7b0eb5368723e4a6d0ccdab07e..088f401945d99a71832d770248e6e721595b3995 100644
--- a/src/species_mod.F90
+++ b/src/species_mod.F90
@@ -29,7 +29,7 @@ MODULE species
REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q2_tau ! factor of the gammaD sum
REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma ! For psi field terms
REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: sqrt_tau_o_sigma ! For Ampere eq
- REAL(xp), ALLOCATABLE, DIMENSION(:),PUBLIC, PROTECTED :: xpdx ! radial pressure gradient
+ REAL(xp), PUBLIC, PROTECTED :: Ptot = 0._dp ! total pressure
!! Accessible routines
PUBLIC :: species_readinputs, species_outputinputs
CONTAINS
@@ -37,7 +37,7 @@ CONTAINS
SUBROUTINE species_readinputs
! Read the input parameters
USE basic, ONLY : lu_in
- USE model, ONLY : Na, nu, ADIAB_E
+ USE model, ONLY : Na, nu, ADIAB_E, MHD_PD
USE prec_const
IMPLICIT NONE
INTEGER :: ia,ib
@@ -47,6 +47,7 @@ CONTAINS
! allocate the arrays of species parameters
CALL species_allocate
! loop over the species namelists in the input file
+ Ptot = 0._xp
DO ia = 1,Na
! default parameters
name_ = 'ions'
@@ -74,13 +75,14 @@ CONTAINS
q2_tau(ia) = (q_**2)/tau_
q_o_sqrt_tau_sigma(ia) = q_/SQRT(tau_)/sigma_
sqrt_tau_o_sigma(ia) = SQRT(tau_)/sigma_
- xpdx(ia) = 0._xp !not implemented yet
+ Ptot = Ptot + tau(ia)*(k_N(ia) + k_T(ia))
! We remove the adiabatic electron flag if electrons are included
SELECT CASE (name_)
CASE ('electrons','e','electron')
ADIAB_E = .FALSE.
END SELECT
ENDDO
+ IF (.NOT. MHD_PD) Ptot = 0._dp
!! Set collision frequency tensor
IF (nu .EQ. 0) THEN
nu_ab = 0
@@ -106,7 +108,6 @@ CONTAINS
! nu_i = nu ! ion-ion collision frequ.
! nu_ee = nu_e ! e-e coll. frequ.
! nu_ie = nu_i ! i-e coll. frequ.
-
END SUBROUTINE species_readinputs
@@ -150,7 +151,6 @@ CONTAINS
ALLOCATE( q2_tau(1:Na))
ALLOCATE(q_o_sqrt_tau_sigma(1:Na))
ALLOCATE( sqrt_tau_o_sigma(1:Na))
- ALLOCATE( xpdx(1:Na))
END SUBROUTINE species_allocate
END MODULE species