From dde677bc010d85121bb5922a7d4f4c50f9bc8856 Mon Sep 17 00:00:00 2001
From: Antoine <antoine.hoffmann@epfl.ch>
Date: Tue, 19 Sep 2023 09:18:53 +0200
Subject: [PATCH] WIP - ExB
---
src/ExB_shear_flow_mod.F90 | 205 +++++++++++++++++++++++--------------
src/auxval.F90 | 2 +-
src/control.F90 | 16 ++-
src/fourier_mod.F90 | 14 +--
src/geometry_mod.F90 | 6 +-
src/grid_mod.F90 | 26 +++--
src/nonlinear_mod.F90 | 6 +-
7 files changed, 167 insertions(+), 108 deletions(-)
diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 73377126..f4804937 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -2,7 +2,7 @@ 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.
- USE prec_const, ONLY: xp, imagu
+ USE prec_const, ONLY: xp, imagu, pi
IMPLICIT NONE
! Variables
@@ -23,9 +23,12 @@ 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 grid, ONLY: Nx, local_nky, total_nky, local_nx, Ny, deltakx, deltaky,&
+ kyarray, kyarray_full
USE geometry, ONLY: Cyq0_x0
+ USE basic, ONLY: dt
IMPLICIT NONE
+ INTEGER :: iky
REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shearing rate and aux var
@@ -42,9 +45,20 @@ CONTAINS
! Setup the ExB shift array
ALLOCATE(sky_ExB(local_nky))
+ ! consider no initial shift (maybe changed if restart)
sky_ExB = 0._xp
- ALLOCATE(sky_ExB_full(total_nky))
+ ! 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)
sky_ExB_full = 0._xp
+ ! 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
! Setup the kx correction array
ALLOCATE(dkx_ExB(local_nky))
@@ -60,27 +74,74 @@ CONTAINS
! Setup nonlinear factor
ALLOCATE( ExB_NL_factor(Nx,local_nky))
- ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx))
+ ALLOCATE(inv_ExB_NL_factor(Ny/2+2,local_nx))
ExB_NL_factor = 1._xp
inv_ExB_NL_factor = 1._xp
END SUBROUTINE Setup_ExB_shear_flow
- ! Update according to the current ExB shear value
- ! -shift the grids
- ! -the spatial operators
- ! -the fields by imposing a shift on kx
+ ! 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, kyarray_full, update_grids, deltaky
+ USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
+ USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
+ USE model, ONLY: LINEARITY
+ USE time_integration, ONLY: c_E
+ IMPLICIT NONE
+ INTEGER, INTENT(IN) :: step_number
+ ! local var
+ INTEGER :: iky
+ REAL(xp):: tnext
+
+ ! tnext = time + c_E(step_number)*dt
+ tnext = time + 0._xp*dt
+ ! do nothing if no ExB
+ IF(ExB) 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
+ 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
+ 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
+ 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, update_grids, &
- total_nkx, deltakx, kx_min, kx_max, kxarray0
+ USE grid, ONLY: local_nky, total_nky, update_grids, &
+ total_nkx, deltakx, kx_min, kx_max, kxarray0, inv_dkx
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
+ INTEGER :: iky, ikx, ikx_s, i_, loopstart, loopend, increment, jump_
IF(ExB) THEN
- ! shift all fields and correct the shift value
+ ! 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
@@ -107,15 +168,12 @@ CONTAINS
ikx = i_ - total_nkx/2 + 1
ENDIF
ikx_s = ikx + jump_ExB(iky)
- ! adjust the shift according
+ ! adjust the shift accordingly
IF (ikx_s .LE. 0) &
ikx_s = ikx_s + total_nkx
IF (ikx_s .GT. total_nkx) &
ikx_s = ikx_s - total_nkx
- ! We test if the shifted modes are still in contained in our resolution
- ! IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
- ! (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
- ! ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
+ ! Then we test if the shifted modes are still in contained in our resolution
IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. kx_max) .AND. &
(kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. kx_min)) THEN
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
@@ -129,73 +187,66 @@ CONTAINS
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.
ENDIF
ENDDO
ENDIF
+ ! 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
- ! update the ExB shear value for the next time step
- SUBROUTINE Update_ExB_shear_flow
- USE basic, ONLY: dt, time
- USE grid, ONLY: local_nky, local_nky_offset, kyarray, inv_dkx, xarray, Nx, Ny, &
- local_nx, local_nx_offset, kyarray_full,&
- ikyarray, inv_ikyarray, deltaky, update_grids
- USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
- USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
+ SUBROUTINE Update_nonlinear_ExB_factors
+ USE grid, ONLY: local_nky, local_nky_offset, xarray, Nx, Ny, local_nx, deltakx,&
+ local_nx_offset, ikyarray, inv_ikyarray, deltaky, update_grids
+ USE basic, ONLY: time, dt
IMPLICIT NONE
- ! local var
INTEGER :: iky, ix
- REAL(xp):: dt_ExB, J_dp, inv_J, x
- ! do nothing if no ExB
- IF(ExB) THEN
- ! reset the ExB shift, jumps and flags
- shiftnow_ExB = .FALSE.
- DO iky = 1,local_Nky
- sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*gamma_E*dt
- 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)
+ REAL(xp):: dt_ExB, J_xp, inv_J, xval, tnext
+ tnext = time + dt
+ 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-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 = (tnext - t0*inv_J*ANINT(J_xp*tnext*inv_t0,xp))
+ DO ix = 1,Nx
+ xval = 2._xp*pi/deltakx*REAL(ix-1,xp)/REAL(Nx,xp)!xarray(ix)
+ ! assemble the ExB nonlin factor
+ ! ExB_NL_factor(ix,iky) = EXP(-imagu*x*gamma_E*J_xp*deltaky*dt_ExB)
+ ExB_NL_factor(ix,iky) = EXP(-imagu*sky_ExB(iky)*xval)
+ ! ExB_NL_factor(ix,iky) = EXP(-imagu*sky_ExB_full(iky+local_nky_offset)*xval)
ENDDO
- ! Update the full skyExB array too
- sky_ExB_full = sky_ExB_full - kyarray_full*gamma_E*dt
- ! Update the ExB nonlinear factor...
- DO iky = 1,local_Nky ! WARNING: Local indices ky loop
- ! for readability
- J_dp = ikyarray(iky+local_nky_offset)
- inv_J = inv_ikyarray(iky+local_nky_offset)
- ! compute dt factor
- dt_ExB = (time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp))
- dkx_ExB(iky) = gamma_E*J_dp*deltaky*dt_ExB
- DO ix = 1,Nx
- x = xarray(ix)
- ! assemble the ExB nonlin factor
- !ExB_NL_factor(ix,iky) = EXP(imagu*x*dkx_ExB(iky))
- ExB_NL_factor(ix,iky) = EXP(imagu*x*sky_ExB(iky)) !GENE does that???
- ENDDO
- ENDDO
- ! ... and the inverse
- DO iky = 1,Ny/2+1 ! WARNING: Global indices ky loop
- ! for readability
- J_dp = ikyarray(iky)
- inv_J = inv_ikyarray(iky)
- ! compute dt factor
- dt_ExB = (time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp))
- DO ix = 1,local_nx
- x = xarray(ix+local_nx_offset)
- ! assemble the inverse ExB nonlin factor
- ! inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*gamma_E*J_dp*deltaky*dt_ExB)
- inv_ExB_NL_factor(iky,ix) = EXP(-imagu*x*sky_ExB_full(iky)) !GENE does that???
- ENDDO
+ ENDDO
+ ! ... and the inverse
+ DO iky = 1,Ny/2+2 ! 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 = (tnext - t0*inv_J*ANINT(J_xp*tnext*inv_t0,xp))
+ DO ix = 1,local_nx
+ xval = 2._xp*pi/deltakx*REAL(ix-1,xp)/REAL(Nx,xp)!xarray(ix+local_nx_offset)
+ ! assemble the inverse ExB nonlin factor
+ ! inv_ExB_NL_factor(iky,ix) = EXP(imagu*x*gamma_E*J_xp*deltaky*dt_ExB)
+ inv_ExB_NL_factor(iky,ix) = EXP(imagu*sky_ExB_full(iky)*xval)
ENDDO
- ENDIF
- ! We update the operators and grids
- ! update the grids
- CALL update_grids(dkx_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
- END SUBROUTINE Update_ExB_shear_flow
+ ENDDO
+ END SUBROUTINE Update_nonlinear_ExB_factors
+
END MODULE ExB_shear_flow
\ No newline at end of file
diff --git a/src/auxval.F90 b/src/auxval.F90
index f5495643..f54b3c38 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -25,7 +25,7 @@ subroutine auxval
CALL speak('=== Set auxiliary values ===')
! Init the grids
CALL init_grids_data(Na,EM,LINEARITY)
- CALL set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
+ CALL set_grids(shear,ExBrate,Npol,LINEARITY,N_HD,EM,Na)
! Allocate memory for global arrays
CALL memory
! Initialize displacement and receive arrays
diff --git a/src/control.F90 b/src/control.F90
index a5379531..a5bf592d 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -10,7 +10,7 @@ SUBROUTINE control
USE initial, ONLY: initialize
USE mpi, ONLY: MPI_COMM_WORLD
USE diagnostics, ONLY: diagnose
- USE ExB_shear_flow, ONLY: Array_shift_ExB_shear_flow, Update_ExB_shear_flow
+ USE ExB_shear_flow, ONLY: Update_ExB_shear_flow
IMPLICIT NONE
REAL(xp) :: t_init_diag_0, t_init_diag_1
INTEGER :: ierr
@@ -67,11 +67,15 @@ SUBROUTINE control
! 2. Main loop
DO
CALL start_chrono(chrono_step) ! Measuring time per step
-
+
! Test if the stopping requirements are met (update nlend)
CALL tesend
IF( nlend ) EXIT ! exit do loop
+ ! Increment steps and csteps (private in basic module)
+ CALL increase_step
+ CALL increase_cstep
+
! Update the ExB shear flow for the next step
! This call includes :
! - the ExB shear value (s(ky)) update for the next time step
@@ -79,15 +83,9 @@ SUBROUTINE control
! - the ExB NL correction factor update (exp(+/- ixkySdts))
! - (optional) the kernel, poisson op. and ampere op update
CALL start_chrono(chrono_ExBs)
- CALL Update_ExB_shear_flow
- ! Shift the arrays if the shear value sky is too high
- CALL Array_shift_ExB_shear_flow
+ CALL Update_ExB_shear_flow(1)
CALL stop_chrono(chrono_ExBs)
- ! Increment steps and csteps (private in basic module)
- CALL increase_step
- CALL increase_cstep
-
! Do a full RK step (e.g. 4 substeps for ERK4)
CALL stepon
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index f50df010..56f1217c 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -266,7 +266,7 @@ END SUBROUTINE fft1D_plans
! module variable (convolution theorem)
SUBROUTINE poisson_bracket_and_sum( ky_array, kx_array, inv_Ny, inv_Nx, AA_y, AA_x,&
local_nky, total_nkx, F_, G_,&
- ExB, ExB_NL_factor, sum_real_)
+ ExB, ExB_NL_factor,sky_ExB,sum_real_)
IMPLICIT NONE
INTEGER, INTENT(IN) :: local_nky,total_nkx
REAL(xp), INTENT(IN) :: inv_Nx, inv_Ny
@@ -274,11 +274,12 @@ END SUBROUTINE fft1D_plans
REAL(xp), DIMENSION(total_nkx), INTENT(IN) :: AA_x
REAL(xp), DIMENSION(local_nky,total_nkx), INTENT(IN) :: kx_array
COMPLEX(c_xp_c), DIMENSION(local_nky,total_nkx), &
- INTENT(IN) :: F_, G_
+ INTENT(IN) :: F_, G_
COMPLEX(xp), DIMENSION(total_nkx,local_nky), &
- INTENT(IN) :: ExB_NL_factor
+ INTENT(IN) :: ExB_NL_factor
LOGICAL, INTENT(IN) :: ExB
- real(c_xp_r), pointer, INTENT(INOUT) :: sum_real_(:,:)
+ REAL(xp),DIMENSION(local_nky), INTENT(IN) :: sky_ExB
+ real(c_xp_r), pointer, INTENT(INOUT) :: sum_real_(:,:)
! local variables
INTEGER :: ikx,iky
COMPLEX(xp), DIMENSION(total_nkx,local_nky) :: ikxF, ikyG, ikyF, ikxG
@@ -330,6 +331,7 @@ END SUBROUTINE fft1D_plans
call fftw_mpi_execute_dft_c2r(planb, ikxG, real_data_g)
#endif
sum_real_ = sum_real_ - real_data_f*real_data_g*inv_Ny*inv_Nx
+
END SUBROUTINE poisson_bracket_and_sum
!******************************************************************************!
@@ -371,7 +373,7 @@ END SUBROUTINE fft1D_plans
SUBROUTINE apply_inv_ExB_NL_factor(fyx,inv_ExB_NL_factor)
IMPLICIT NONE
real(c_xp_r), pointer, INTENT(INOUT) :: fyx(:,:)
- COMPLEX(xp), DIMENSION(NY_halved,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
+ COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_), INTENT(IN) :: inv_ExB_NL_factor
! local variables
REAL(xp), DIMENSION(2*NY_halved,local_nx_) :: tmp_yx_1, tmp_yx_2
COMPLEX(xp), DIMENSION(NY_halved+1,local_nx_) :: tmp_kyx
@@ -385,7 +387,7 @@ END SUBROUTINE fft1D_plans
! Fourier real to complex on the second buffer (first buffer is now unusable)
CALL FFT_yx_to_kyx(tmp_yx_1,tmp_kyx)
! Treat the result with the ExB NL factor
- DO iky = 1,NY_halved
+ DO iky = 1,NY_halved+1
DO ix = 1,local_nx_
tmp_kyx(iky,ix) = tmp_kyx(iky,ix)*inv_ExB_NL_factor(iky,ix)
ENDDO
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 91b83db3..73b4e589 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -104,7 +104,7 @@ CONTAINS
USE miller, ONLY: set_miller_parameters, get_miller
USE circular, ONLY: get_circ
USE calculus, ONLY: simpson_rule_z
- USE model, ONLY: beta, LINEARITY, N_HD
+ USE model, ONLY: beta, ExBrate, LINEARITY, N_HD
USE species, ONLY: Ptot
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
@@ -127,6 +127,7 @@ CONTAINS
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
+ C_y = 1._xp
Cyq0_x0 = 1._xp
CASE('z-pinch','Z-pinch','zpinch','Zpinch')
CALL speak('Z-pinch geometry')
@@ -136,6 +137,7 @@ CONTAINS
q0 = 0._xp
eps = 0._xp
kappa = 1._xp
+ C_y = 1._xp
Cyq0_x0 = 1._xp
CASE('miller','Miller')
CALL speak('Miller geometry')
@@ -165,7 +167,7 @@ CONTAINS
! 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)
+ CALL set_kxgrid(shear,ExBrate,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}
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index a17d89df..cec77fa8 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -292,8 +292,8 @@ CONTAINS
! The other grids are simply
! |_|_|_|_|
! 1 2 3 4
- SUBROUTINE set_grids(shear,Npol,LINEARITY,N_HD,EM,Na)
- REAL(xp), INTENT(IN) :: shear, Npol
+ SUBROUTINE set_grids(shear,ExBrate,Npol,LINEARITY,N_HD,EM,Na)
+ REAL(xp), INTENT(IN) :: shear, ExBrate, Npol
CHARACTER(len=*), INTENT(IN) :: LINEARITY
INTEGER, INTENT(IN) :: N_HD
LOGICAL, INTENT(IN) :: EM
@@ -302,7 +302,7 @@ CONTAINS
CALL set_pgrid(EM)
CALL set_jgrid
CALL set_kygrid(LINEARITY,N_HD)
- CALL set_kxgrid(shear,Npol,1._xp,LINEARITY,N_HD) ! this will be redone after geometry if Cyq0_x0 .NE. 1
+ CALL set_kxgrid(shear,ExBrate,Npol,1._xp,LINEARITY,N_HD) ! this will be redone after geometry if Cyq0_x0 .NE. 1
CALL set_xgrid
CALL set_zgrid (Npol)
END SUBROUTINE set_grids
@@ -486,10 +486,10 @@ CONTAINS
ENDIF
END SUBROUTINE set_kygrid
- SUBROUTINE set_kxgrid(shear,Npol,Cyq0_x0,LINEARITY,N_HD)
- USE prec_const
+ SUBROUTINE set_kxgrid(shear,ExBrate,Npol,Cyq0_x0,LINEARITY,N_HD)
+ USE prec_const, ONLY: xp, pi
IMPLICIT NONE
- REAL(xp), INTENT(IN) :: shear, Npol, Cyq0_x0
+ REAL(xp), INTENT(IN) :: shear, Npol, Cyq0_x0, ExBrate
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
INTEGER :: ikx, iky
@@ -538,7 +538,13 @@ CONTAINS
ERROR STOP "Gyacomo is safer with an even Kx number"
ENDIF
! Orszag 2/3 filter
- two_third_kxmax = 2._xp/3._xp*(kx_max-deltakx);
+ ! We remove one point more if ExB is on since the moving grid
+ ! can go up to kx=+-(kx_max+deltakx/2)
+ IF (ABS(ExBrate) .GT. 0) THEN
+ two_third_kxmax = 2._xp/3._xp*(kx_max-deltakx)
+ ELSE
+ two_third_kxmax = 2._xp/3._xp*kx_max
+ ENDIF
! Antialiasing filter
DO iky = 1,local_nky
DO ikx = 1,local_nkx
@@ -663,16 +669,16 @@ CONTAINS
two_third_kpmax = 2._xp/3._xp * MAXVAL(kparray)
END SUBROUTINE
- SUBROUTINE update_grids (dkx_ExB,gxx,gxy,gyy,inv_hatB2)
+ SUBROUTINE update_grids (sky_ExB,gxx,gxy,gyy,inv_hatB2)
IMPLICIT NONE
- REAL(xp), DIMENSION(local_nky),INTENT(IN) :: dkx_ExB ! ExB correction dkx = gamma_E ky dtshift
+ REAL(xp), DIMENSION(local_nky),INTENT(IN) :: sky_ExB ! ExB correction dkx = gamma_E ky dtshift
REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
INTEGER :: eo,iz,iky,ikx
REAL(xp) :: kx, ky
! Update the kx grid
DO ikx = 1,total_Nkx
DO iky = 1,local_nky
- kxarray(iky,ikx) = kxarray0(ikx) - dkx_ExB(iky)
+ kxarray(iky,ikx) = kxarray0(ikx) - sky_ExB(iky)
ENDDO
ENDDO
! Update the kperp grid
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 946304ac..2440078e 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -14,7 +14,7 @@ MODULE nonlinear
USE prec_const, ONLY : xp
USE species, ONLY : sqrt_tau_o_sigma
USE time_integration, ONLY : updatetlevel
- USE ExB_shear_flow, ONLY : ExB_NL_factor, inv_ExB_NL_factor, ExB, dkx_ExB
+ USE ExB_shear_flow, ONLY : ExB_NL_factor, inv_ExB_NL_factor, ExB, sky_ExB
use, intrinsic :: iso_c_binding
IMPLICIT NONE
@@ -100,7 +100,7 @@ SUBROUTINE compute_nonlinear
! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky,total_nkx,F_cmpx,G_cmpx,&
- ExB, ExB_NL_factor, bracket_sum_r)
+ ExB, ExB_NL_factor, sky_ExB, bracket_sum_r)
!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
IF(EM) THEN
! First convolution terms
@@ -116,7 +116,7 @@ SUBROUTINE compute_nonlinear
! this function adds its result to bracket_sum_r
CALL poisson_bracket_and_sum( kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,&
local_nky,total_nkx,F_cmpx,G_cmpx,&
- ExB, ExB_NL_factor, bracket_sum_r)
+ ExB, ExB_NL_factor, sky_ExB, bracket_sum_r)
ENDIF
ENDDO n
! Apply the ExB shearing rate factor before going back to k-space
--
GitLab