Newer
Older
MODULE ExB_shear_flow
! This module contains the necessary tools to implement ExB shearing flow effects.
! The algorithm is taken from the presentation of Hammett et al. 2006 (APS) and
! it the one used in GS2.
IMPLICIT NONE
! Variables
REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E
REAL(xp), PUBLIC, PROTECTED :: t0, inv_t0 = 0._xp ! charact. shear time
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky)
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB_full ! full ky version
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jump_ExB ! jump to do to shift the kx grids
LOGICAL, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: ExB_NL_factor! factor for nonlinear term
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: inv_ExB_NL_factor
LOGICAL, PUBLIC, PROTECTED :: ExB = .false. ! presence of ExB background shearing rate
REAL(xp), PUBLIC, PROTECTED :: maximal_kx, minimal_kx ! max and min kx evolved for array shifting
! Routines
PUBLIC :: Setup_ExB_shear_flow, Array_shift_ExB_shear_flow, Update_ExB_shear_flow
CONTAINS
! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow(ExBrate)
USE grid, ONLY: Nx, local_nky, total_nky, local_nx, Ny, deltakx, deltaky,&
USE geometry, ONLY: Cyq0_x0, C_y
USE model, ONLY: LINEARITY
IMPLICIT NONE
! Setup the ExB shearing rate and aux var
! In GENE, there is a minus sign here...
gamma_E = ExBrate*C_y*abs(Cyq0_x0/C_y)
IF(abs(gamma_E) .GT. EPSILON(gamma_E)) THEN
CALL speak('-ExB background flow detected-',2)
ExB = .TRUE.
t0 = deltakx/deltaky/gamma_E
inv_t0 = 1._xp/t0
ELSE ! avoid 1/0 division (t0 is killed anyway in this case)
ExB = .FALSE.
t0 = 0._xp
inv_t0 = 0._xp
ENDIF
! Setup the ExB shift array
ALLOCATE(sky_ExB(local_nky))
! consider no initial shift (maybe changed if restart)
sky_ExB = 0._xp
! Midpoint init
DO iky = 1,local_nky
sky_ExB(iky) = sky_ExB(iky) !+ 0.5_xp*kyarray(iky)*gamma_E*dt
ENDDO
ALLOCATE(sky_ExB_full(total_nky+1))
! consider no initial shift (maybe changed if restart)
! Midpoint init
DO iky = 1,total_nky+1
sky_ExB_full(iky) = sky_ExB_full(iky) !+ 0.5_xp*REAL(iky-1,xp)*deltaky*gamma_E*dt
ENDDO
ALLOCATE(jump_ExB(local_nky))
jump_ExB = 0
! Setup the shifting flag array
ALLOCATE(shiftnow_ExB(local_nky))
shiftnow_ExB = .FALSE.
! Setup nonlinear factor (McMillan 2019)
ALLOCATE( ExB_NL_factor(Nx,local_nky))
ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx))
ExB_NL_factor = 1._xp
inv_ExB_NL_factor = 1._xp
! Setup maximal evolved modes for the shift conditions
SELECT CASE (LINEARITY)
CASE('linear') ! If linear we just use the max kx
maximal_kx = kx_max
minimal_kx = kx_min
CASE('nonlinear') ! If NL, 2/3 rule
maximal_kx = 2._xp/3._xp*kx_max-deltakx
minimal_kx = 2._xp/3._xp*kx_min+deltakx
END SELECT
END SUBROUTINE Setup_ExB_shear_flow
! update the ExB shear value for the next time step
SUBROUTINE Update_ExB_shear_flow(step_number)
USE basic, ONLY: dt!,time
USE grid, ONLY: local_nky, total_nky, kyarray, inv_dkx, update_grids, deltaky!,kyarray_full
USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
USE time_integration, ONLY: c_E!, ntimelevel
IMPLICIT NONE
INTEGER, INTENT(IN) :: step_number
! local var
INTEGER :: iky
! do nothing if no ExB
IF(ExB) THEN
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
IF (step_number .LT. 0) THEN ! not updated each substeps (call from control)
dt_sub = dt
ELSEIF(step_number .GT. 1) THEN ! updated each substeps
dt_sub = (c_E(step_number)-c_E(step_number-1))*dt
ELSE ! first step is at t=0
dt_sub = 0._xp
ENDIF
! Do nothing if dt is 0
IF(dt_sub .GT. epsilon(dt_sub)) THEN
! Update new shear value
DO iky = 1,local_nky
!! This must be done incrementely to be able to pull it back
! when a grid shift occurs
sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*gamma_E*dt_sub
jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx)
! If the jump is 1 or more for a given ky, we flag the index
! in shiftnow_ExB and will use it in Shift_fields to avoid
! zero-shiftings that may be majoritary.
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
ENDDO
! Update the full skyExB array too
DO iky = 1,total_nky+1
sky_ExB_full(iky) = sky_ExB_full(iky) - REAL(iky-1,xp)*deltaky*gamma_E*dt_sub
ENDDO
! Shift the arrays if the shear value sky is too high
CALL Array_shift_ExB_shear_flow
! We update the operators and grids
! update the grids
CALL update_grids(sky_ExB,gxx,gxy,gyy,inv_hatB2)
! update the EM op., the kernels and the curvature op.
CALL evaluate_kernels
CALL evaluate_EM_op
CALL evaluate_magn_curv
! update the ExB nonlinear factor...
! IF(LINEARITY .EQ. 'nonlinear') &
! CALL Update_nonlinear_ExB_factors(dt_sub)
ENDIF
END SUBROUTINE Update_ExB_shear_flow
! According to the current ExB shear value we update
! the fields by imposing a shift on kx
SUBROUTINE Array_shift_ExB_shear_flow
USE grid, ONLY: local_nky, total_nky, update_grids, &
total_nkx, deltakx, kxarray0, inv_dkx!,kx_min, kx_max
USE prec_const, ONLY: PI
USE fields, ONLY: moments, phi, psi
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
IMPLICIT NONE
! local var
INTEGER :: iky, ikx, ikx_s, i_, loopstart, loopend, increment, jump_
! shift all local fields and correct the local shift value
DO iky = 1,local_Nky
IF(shiftnow_ExB(iky)) THEN
! We shift the array from left to right or right to left according to the jump
! This avoids to make copy
IF(jump_ExB(iky) .LT. 0) THEN
loopstart = 1
loopend = total_nkx
increment = 1
ELSE
loopstart = total_nkx
loopend = 1
increment = -1
!loop to go through the array in a monotonic kx order
! Recall: the kx array is organized as
! 6 7 8 1 2 3 4 5 (indices ikx)
! -3 -2 -1 0 1 2 3 4 (values in dkx)
! to monotonically travel across the kx array the indices write
! 67812345 for positive shift or 54321678 for negative shift
DO i_ = loopstart, loopend, increment
! We shift our index since kx is stored in a [0 kmax]U[kmin -dk] fashion
IF (i_ .LT. total_nkx/2) THEN ! go to the negative kx region
ikx = i_ + total_nkx/2 + 1
ELSE ! positive
ikx = i_ - total_nkx/2 + 1
ENDIF
ikx_s = ikx - jump_ExB(iky)
IF (ikx_s .LE. 0) &
ikx_s = ikx_s + total_nkx
IF (ikx_s .GT. total_nkx) &
ikx_s = ikx_s - total_nkx
! Then we test if the shifted modes are still contained in our resolution
IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. maximal_kx) .AND. &
(kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. minimal_kx)) THEN
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
phi(iky,ikx,:) = phi(iky,ikx_s,:)
psi(iky,ikx,:) = psi(iky,ikx_s,:)
ELSE ! if it is not, it is lost (~dissipation for high modes)
moments(:,:,:,iky,ikx,:,:) = 0._xp
phi(iky,ikx,:) = 0._xp
psi(iky,ikx,:) = 0._xp
ENDIF
ENDDO
! correct the shift value s(ky) for this row
sky_ExB(iky) = sky_ExB(iky) - jump_ExB(iky)*deltakx
! reset the flag
shiftnow_ExB(iky) = .FALSE.
! Check the global shift values
DO iky = 1,total_nky+1
jump_ = NINT(sky_ExB_full(iky)*inv_dkx)
IF (ABS(jump_) .GT. 0) &
sky_ExB_full(iky) = sky_ExB_full(iky) - jump_*deltakx
ENDDO
END SUBROUTINE Array_shift_ExB_shear_flow
SUBROUTINE Update_nonlinear_ExB_factors(dt_sub)
USE grid, ONLY: local_nky, local_nky_offset, Nx, Ny, local_nx, deltakx,&
local_nx_offset, deltaky, update_grids!,xarray, ikyarray, inv_ikyarray
USE basic, ONLY: time!, dt
! USE time_integration, ONLY: c_E
REAL(xp), INTENT(IN) :: dt_sub
INTEGER :: iky, ix
REAL(xp):: dt_ExB, J_xp, inv_J, &
I_xp, deltax, Lx, xval, tnow, kxExB, dkx_ExB, v_ExB
! aux val
tnow = time + dt_sub
Lx = 2._xp*pi/deltakx
deltax = Lx/REAL(Nx,xp)
DO iky = 1,local_nky ! WARNING: Local indices ky loop
! for readability
! J_xp = ikyarray(iky+local_nky_offset)
! inv_J = inv_ikyarray(iky+local_nky_offset)
J_xp = REAL(iky+local_nky_offset-1,xp)
IF(J_xp .GT. 0._xp) THEN
inv_J = 1._xp/J_xp
ELSE
inv_J = 0._xp
ENDIF
! compute dt factor
dt_ExB = (tnow - t0*inv_J*ANINT(J_xp*tnow*inv_t0,xp))
v_ExB = -gamma_E*J_xp*deltaky
dkx_ExB= v_ExB*tnow - ANINT(v_ExB*tnow/deltakx,xp)*deltakx
I_xp = REAL(ix-1,xp)
xval = I_xp*deltax
kxExB = gamma_E*J_xp*deltaky*dt_ExB
! kxExB = sky_ExB(iky) ! GENE
kxExB = dkx_ExB ! home made
ExB_NL_factor(ix,iky) = EXP(-imagu*kxExB*xval)
! ExB_NL_factor(ix,iky) = 1._xp + (imagu*kxExB*xval)
DO iky = 1,Ny/2+1 ! WARNING: Global indices ky loop
! for readability
J_xp = REAL(iky-1,xp)
IF(J_xp .GT. 0._xp) THEN
inv_J = 1._xp/J_xp
ELSE
inv_J = 0._xp
ENDIF
! compute dt factor
dt_ExB = tnow - t0*inv_J*ANINT(J_xp*tnow*inv_t0,xp)
v_ExB = -gamma_E*J_xp*deltaky
dkx_ExB= v_ExB*tnow - ANINT(v_ExB*tnow/deltakx,xp)*deltakx
I_xp = REAL(ix+local_nx_offset-1,xp)
xval = I_xp*deltax
! kxExB = gamma_E*J_xp*deltaky*dt_ExB
! kxExB = sky_ExB_full(iky) ! GENE
kxExB = dkx_ExB ! home made
inv_ExB_NL_factor(iky,ix) = EXP(imagu*kxExB*xval)
! inv_ExB_NL_factor(iky,ix) = 1._xp + (imagu*kxExB*xval)
END SUBROUTINE Update_nonlinear_ExB_factors
END MODULE ExB_shear_flow