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

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
parent 5348e128
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
+ 101
− 63
View file @ d5560ced
...@@ -10,12 +10,12 @@ MODULE ExB_shear_flow ...@@ -10,12 +10,12 @@ MODULE ExB_shear_flow
REAL(xp), PUBLIC, PROTECTED :: t0, inv_t0 = 0._xp ! charact. shear time 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 ! 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 :: 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 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 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 :: ExB_NL_factor! factor for nonlinear term
COMPLEX(xp),DIMENSION(:,:), ALLOCATABLE, PUBLIC, PROTECTED :: inv_ExB_NL_factor 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 ! Routines
PUBLIC :: Setup_ExB_shear_flow, Array_shift_ExB_shear_flow, Update_ExB_shear_flow PUBLIC :: Setup_ExB_shear_flow, Array_shift_ExB_shear_flow, Update_ExB_shear_flow
...@@ -24,15 +24,17 @@ CONTAINS ...@@ -24,15 +24,17 @@ CONTAINS
! Setup the variables for the ExB shear ! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow(ExBrate) 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 kyarray, kyarray_full, kx_max, kx_min
USE geometry, ONLY: Cyq0_x0 USE geometry, ONLY: Cyq0_x0, C_y
USE basic, ONLY: dt USE basic, ONLY: dt
USE model, ONLY: LINEARITY
IMPLICIT NONE IMPLICIT NONE
INTEGER :: iky INTEGER :: iky
REAL(xp), INTENT(IN) :: ExBrate REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shearing rate and aux var ! 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 IF(abs(gamma_E) .GT. EPSILON(gamma_E)) THEN
ExB = .TRUE. ExB = .TRUE.
t0 = deltakx/deltaky/gamma_E t0 = deltakx/deltaky/gamma_E
...@@ -60,10 +62,6 @@ CONTAINS ...@@ -60,10 +62,6 @@ CONTAINS
sky_ExB_full(iky) = sky_ExB_full(iky) !+ 0.5_xp*REAL(iky-1,xp)*deltaky*gamma_E*dt sky_ExB_full(iky) = sky_ExB_full(iky) !+ 0.5_xp*REAL(iky-1,xp)*deltaky*gamma_E*dt
ENDDO ENDDO
! Setup the kx correction array
ALLOCATE(dkx_ExB(local_nky))
dkx_ExB = 0._xp
! Setup the jump array ! Setup the jump array
ALLOCATE(jump_ExB(local_nky)) ALLOCATE(jump_ExB(local_nky))
jump_ExB = 0 jump_ExB = 0
...@@ -78,6 +76,16 @@ CONTAINS ...@@ -78,6 +76,16 @@ CONTAINS
ExB_NL_factor = 1._xp ExB_NL_factor = 1._xp
inv_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 END SUBROUTINE Setup_ExB_shear_flow
! update the ExB shear value for the next time step ! update the ExB shear value for the next time step
...@@ -87,42 +95,52 @@ CONTAINS ...@@ -87,42 +95,52 @@ CONTAINS
USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv USE geometry, ONLY: gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
USE numerics, ONLY: evaluate_EM_op, evaluate_kernels USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
USE model, ONLY: LINEARITY USE model, ONLY: LINEARITY
USE time_integration, ONLY: c_E USE time_integration, ONLY: c_E, ntimelevel
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(IN) :: step_number INTEGER, INTENT(IN) :: step_number
! local var ! local var
REAL(xp):: dt_sub
INTEGER :: iky INTEGER :: iky
! do nothing if no ExB ! do nothing if no ExB
IF(ExB) THEN IF(ExB) THEN
! Update new shear value IF (step_number .LT. 0) THEN ! not updated each substeps (call from control)
DO iky = 1,local_nky dt_sub = dt
!! This must be done incrementely to be able to pull it back ELSEIF(step_number .GT. 1) THEN ! updated each substeps
! when a grid shift occurs dt_sub = (c_E(step_number)-c_E(step_number-1))*dt
sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*gamma_E*dt ELSE ! first step is at t=0
jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx) dt_sub = 0._xp
! If the jump is 1 or more for a given ky, we flag the index ENDIF
! in shiftnow_ExB and will use it in Shift_fields to avoid ! Do nothing if dt is 0
! zero-shiftings that may be majoritary. IF(dt_sub .GT. epsilon(dt_sub)) THEN
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0) ! Update new shear value
ENDDO DO iky = 1,local_nky
! Update the full skyExB array too !! This must be done incrementely to be able to pull it back
DO iky = 1,total_nky+1 ! when a grid shift occurs
sky_ExB_full(iky) = sky_ExB_full(iky) - REAL(iky-1,xp)*deltaky*gamma_E*dt sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*gamma_E*dt_sub
ENDDO jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx)
! Shift the arrays if the shear value sky is too high ! If the jump is 1 or more for a given ky, we flag the index
CALL Array_shift_ExB_shear_flow ! in shiftnow_ExB and will use it in Shift_fields to avoid
! zero-shiftings that may be majoritary.
! We update the operators and grids shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
! update the grids ENDDO
CALL update_grids(sky_ExB,gxx,gxy,gyy,inv_hatB2) ! Update the full skyExB array too
! update the EM op., the kernels and the curvature op. DO iky = 1,total_nky+1
CALL evaluate_kernels sky_ExB_full(iky) = sky_ExB_full(iky) - REAL(iky-1,xp)*deltaky*gamma_E*dt_sub
CALL evaluate_EM_op ENDDO
CALL evaluate_magn_curv ! Shift the arrays if the shear value sky is too high
! update the ExB nonlinear factor... CALL Array_shift_ExB_shear_flow
IF(LINEARITY .EQ. 'nonlinear') &
CALL update_nonlinear_ExB_factors(step_number) ! 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 ENDIF
END SUBROUTINE Update_ExB_shear_flow END SUBROUTINE Update_ExB_shear_flow
...@@ -143,7 +161,7 @@ CONTAINS ...@@ -143,7 +161,7 @@ CONTAINS
IF(shiftnow_ExB(iky)) THEN IF(shiftnow_ExB(iky)) THEN
! We shift the array from left to right or right to left according to the jump ! We shift the array from left to right or right to left according to the jump
! This avoids to make copy ! This avoids to make copy
IF(jump_ExB(iky) .GT. 0) THEN IF(jump_ExB(iky) .LT. 0) THEN
loopstart = 1 loopstart = 1
loopend = total_nkx loopend = total_nkx
increment = 1 increment = 1
...@@ -156,23 +174,24 @@ CONTAINS ...@@ -156,23 +174,24 @@ CONTAINS
! Recall: the kx array is organized as ! Recall: the kx array is organized as
! 6 7 8 1 2 3 4 5 (indices ikx) ! 6 7 8 1 2 3 4 5 (indices ikx)
! -3 -2 -1 0 1 2 3 4 (values in dkx) ! -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 ! to monotonically travel across the kx array the indices write
! 67812345 or 54321678 ! 67812345 for positive shift or 54321678 for negative shift
DO i_ = loopstart, loopend, increment 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 IF (i_ .LT. total_nkx/2) THEN ! go to the negative kx region
ikx = i_ + total_nkx/2 + 1 ikx = i_ + total_nkx/2 + 1
ELSE ! positive ELSE ! positive
ikx = i_ - total_nkx/2 + 1 ikx = i_ - total_nkx/2 + 1
ENDIF ENDIF
ikx_s = ikx + jump_ExB(iky) ikx_s = ikx - jump_ExB(iky)
! adjust the shift accordingly ! adjust the shift accordingly
IF (ikx_s .LE. 0) & IF (ikx_s .LE. 0) &
ikx_s = ikx_s + total_nkx ikx_s = ikx_s + total_nkx
IF (ikx_s .GT. total_nkx) & IF (ikx_s .GT. total_nkx) &
ikx_s = ikx_s - total_nkx ikx_s = ikx_s - total_nkx
! Then we test if the shifted modes are still in contained in our resolution ! Then we test if the shifted modes are still contained in our resolution
IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. kx_max) .AND. & IF ( (kxarray0(ikx)+jump_ExB(iky)*deltakx .LE. maximal_kx) .AND. &
(kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. kx_min)) THEN (kxarray0(ikx)+jump_ExB(iky)*deltakx .GE. minimal_kx)) THEN
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:) moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
phi(iky,ikx,:) = phi(iky,ikx_s,:) phi(iky,ikx,:) = phi(iky,ikx_s,:)
psi(iky,ikx,:) = psi(iky,ikx_s,:) psi(iky,ikx,:) = psi(iky,ikx_s,:)
...@@ -198,23 +217,26 @@ CONTAINS ...@@ -198,23 +217,26 @@ CONTAINS
END SUBROUTINE Array_shift_ExB_shear_flow 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,& USE grid, ONLY: local_nky, local_nky_offset, xarray, Nx, Ny, local_nx, deltakx,&
local_nx_offset, ikyarray, inv_ikyarray, deltaky, update_grids local_nx_offset, ikyarray, inv_ikyarray, deltaky, update_grids
USE basic, ONLY: time, dt USE basic, ONLY: time, dt
USE time_integration, ONLY: c_E USE time_integration, ONLY: c_E
IMPLICIT NONE IMPLICIT NONE
INTEGER, INTENT(IN) :: step_number REAL(xp), INTENT(IN) :: dt_sub
INTEGER :: iky, ix INTEGER :: iky, ix
REAL(xp):: dt_ExB, J_xp, inv_J, xval, tnow REAL(xp):: dt_ExB, J_xp, inv_J, &
I_xp, deltax, Lx, xval, tnow, kxExB, dkx_ExB, v_ExB
tnow = time + c_E(step_number)*dt ! 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 DO iky = 1,local_nky ! WARNING: Local indices ky loop
! for readability ! for readability
! J_xp = ikyarray(iky+local_nky_offset) ! J_xp = ikyarray(iky+local_nky_offset)
! inv_J = inv_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 IF(J_xp .GT. 0._xp) THEN
inv_J = 1._xp/J_xp inv_J = 1._xp/J_xp
ELSE ELSE
...@@ -222,12 +244,16 @@ CONTAINS ...@@ -222,12 +244,16 @@ CONTAINS
ENDIF ENDIF
! compute dt factor ! 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,Nx DO ix = 1,Nx
xval = 2._xp*pi/deltakx*REAL(ix-1,xp)/REAL(Nx,xp)!xarray(ix) I_xp = REAL(ix-1,xp)
! assemble the ExB nonlin factor xval = I_xp*deltax
ExB_NL_factor(ix,iky) = EXP(-imagu*xval*gamma_E*J_xp*deltaky*dt_ExB) kxExB = gamma_E*J_xp*deltaky*dt_ExB
! ExB_NL_factor(ix,iky) = EXP(-imagu*sky_ExB(iky)*xval) ! kxExB = sky_ExB(iky) ! GENE
! ExB_NL_factor(ix,iky) = EXP(-imagu*sky_ExB_full(iky+local_nky_offset)*xval) 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
ENDDO ENDDO
! ... and the inverse ! ... and the inverse
...@@ -240,14 +266,26 @@ CONTAINS ...@@ -240,14 +266,26 @@ CONTAINS
inv_J = 0._xp inv_J = 0._xp
ENDIF ENDIF
! compute dt factor ! 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 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 ! 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*kxExB*xval)
! inv_ExB_NL_factor(iky,ix) = EXP(imagu*sky_ExB_full(iky)*xval) ! inv_ExB_NL_factor(iky,ix) = 1._xp + (imagu*kxExB*xval)
ENDDO ENDDO
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 END MODULE ExB_shear_flow
\ No newline at end of file
src/control.F90
+ 1
− 1
View file @ d5560ced
...@@ -83,7 +83,7 @@ SUBROUTINE control ...@@ -83,7 +83,7 @@ SUBROUTINE control
! - the ExB NL correction factor update (exp(+/- ixkySdts)) ! - the ExB NL correction factor update (exp(+/- ixkySdts))
! - (optional) the kernel, poisson op. and ampere op update ! - (optional) the kernel, poisson op. and ampere op update
CALL start_chrono(chrono_ExBs) CALL start_chrono(chrono_ExBs)
CALL Update_ExB_shear_flow(1) CALL Update_ExB_shear_flow(-1)
CALL stop_chrono(chrono_ExBs) CALL stop_chrono(chrono_ExBs)
! Do a full RK step (e.g. 4 substeps for ERK4) ! Do a full RK step (e.g. 4 substeps for ERK4)
......
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