Skip to content
Snippets Groups Projects
Commit dde677bc authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

WIP - ExB

parent ff54cc94
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ExB_shear_flow_mod.F90
+ 128
− 77
View file @ dde677bc
......@@ -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
src/auxval.F90
+ 1
− 1
View file @ dde677bc
......@@ -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
......
src/control.F90
+ 7
− 9
View file @ dde677bc
......@@ -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
......
src/fourier_mod.F90
+ 8
− 6
View file @ dde677bc
......@@ -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
......
src/geometry_mod.F90
+ 4
− 2
View file @ dde677bc
......@@ -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}
......
src/grid_mod.F90
+ 16
− 10
View file @ dde677bc
......@@ -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
......
src/nonlinear_mod.F90
+ 3
− 3
View file @ dde677bc
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment