From d5560ced671a0330ee69daba78cff0c6e9d0f3e2 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 26 Sep 2023 16:17:13 +0200
Subject: [PATCH] Update ExB correction (not working) -Update_ExB_shear_flow
asks a step num -this step num can be set to -1 when used in control -the
sign have been checked -coherent with Hamett et al. -the NL factor comp. is
not optmized but fully detailed -still has issue to verify McMillan et al.
test
---
src/ExB_shear_flow_mod.F90 | 164 +++++++++++++++++++++++--------------
src/control.F90 | 2 +-
2 files changed, 102 insertions(+), 64 deletions(-)
diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index b165cb73..97e1fabf 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -10,12 +10,12 @@ MODULE ExB_shear_flow
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
- REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: dkx_ExB ! correction to obtain the exact kx mode
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
+ 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
@@ -24,15 +24,17 @@ 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,&
- kyarray, kyarray_full
- USE geometry, ONLY: Cyq0_x0
+ kyarray, kyarray_full, kx_max, kx_min
+ USE geometry, ONLY: Cyq0_x0, C_y
USE basic, ONLY: dt
+ USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: iky
REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shearing rate and aux var
- gamma_E = -ExBrate*Cyq0_x0
+ ! 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
ExB = .TRUE.
t0 = deltakx/deltaky/gamma_E
@@ -60,10 +62,6 @@ CONTAINS
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))
- dkx_ExB = 0._xp
-
! Setup the jump array
ALLOCATE(jump_ExB(local_nky))
jump_ExB = 0
@@ -78,6 +76,16 @@ CONTAINS
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
@@ -87,42 +95,52 @@ CONTAINS
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
+ USE time_integration, ONLY: c_E, ntimelevel
IMPLICIT NONE
INTEGER, INTENT(IN) :: step_number
! local var
+ REAL(xp):: dt_sub
INTEGER :: iky
-
! 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(step_number)
+ 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
ENDIF
END SUBROUTINE Update_ExB_shear_flow
@@ -143,7 +161,7 @@ CONTAINS
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) .GT. 0) THEN
+ IF(jump_ExB(iky) .LT. 0) THEN
loopstart = 1
loopend = total_nkx
increment = 1
@@ -156,23 +174,24 @@ CONTAINS
! 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)
- ! so to go along the array in a monotonic way one must travel as
- ! 67812345 or 54321678
+ ! 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)
+ ikx_s = ikx - jump_ExB(iky)
! 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
- ! 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
+ ! 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,:)
@@ -198,23 +217,26 @@ CONTAINS
END SUBROUTINE Array_shift_ExB_shear_flow
- SUBROUTINE Update_nonlinear_ExB_factors(step_number)
+ SUBROUTINE Update_nonlinear_ExB_factors(dt_sub)
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
USE time_integration, ONLY: c_E
IMPLICIT NONE
- INTEGER, INTENT(IN) :: step_number
+ REAL(xp), INTENT(IN) :: dt_sub
INTEGER :: iky, ix
- REAL(xp):: dt_ExB, J_xp, inv_J, xval, tnow
-
- tnow = time + c_E(step_number)*dt
+ 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-1,xp)
+ J_xp = REAL(iky+local_nky_offset-1,xp)
IF(J_xp .GT. 0._xp) THEN
inv_J = 1._xp/J_xp
ELSE
@@ -222,12 +244,16 @@ CONTAINS
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
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*xval*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)
+ 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)
ENDDO
ENDDO
! ... and the inverse
@@ -240,14 +266,26 @@ CONTAINS
inv_J = 0._xp
ENDIF
! compute dt factor
- dt_ExB = (tnow - t0*inv_J*ANINT(J_xp*tnow*inv_t0,xp))
+ 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
DO ix = 1,local_nx
- xval = 2._xp*pi/deltakx*REAL(ix-1,xp)/REAL(Nx,xp)!xarray(ix+local_nx_offset)
+ 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
! assemble the inverse ExB nonlin factor
- inv_ExB_NL_factor(iky,ix) = EXP(imagu*gamma_E*J_xp*deltaky*dt_ExB*xval)
- ! inv_ExB_NL_factor(iky,ix) = EXP(imagu*sky_ExB_full(iky)*xval)
+ inv_ExB_NL_factor(iky,ix) = EXP(imagu*kxExB*xval)
+ ! inv_ExB_NL_factor(iky,ix) = 1._xp + (imagu*kxExB*xval)
ENDDO
ENDDO
- END SUBROUTINE Update_nonlinear_ExB_factors
-
+ ! test
+ DO iky = 1,Ny/2+1
+ DO ix = 1,Nx
+ write(42,*) REAL(inv_ExB_NL_factor(iky,ix) * ExB_NL_factor(ix,iky))
+ ENDDO
+ ENDDO
+ stop
+ END SUBROUTINE Update_nonlinear_ExB_factors
END MODULE ExB_shear_flow
\ No newline at end of file
diff --git a/src/control.F90 b/src/control.F90
index a5bf592d..0095c2eb 100644
--- a/src/control.F90
+++ b/src/control.F90
@@ -83,7 +83,7 @@ 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(1)
+ CALL Update_ExB_shear_flow(-1)
CALL stop_chrono(chrono_ExBs)
! Do a full RK step (e.g. 4 substeps for ERK4)
--
GitLab