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

move the ExB flag

parent a16d890f
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
+ 113
− 108
View file @ 4b80c283
...@@ -13,24 +13,39 @@ MODULE ExB_shear_flow ...@@ -13,24 +13,39 @@ MODULE ExB_shear_flow
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
! Routines ! Routines
PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
CONTAINS CONTAINS
! Setup the variables for the ExB shear ! Setup the variables for the ExB shear
SUBROUTINE Setup_ExB_shear_flow SUBROUTINE Setup_ExB_shear_flow(ExBrate)
USE grid, ONLY : Nx, local_nky, local_nx, Ny, deltakx, deltaky USE grid, ONLY : Nx, local_nky, local_nx, Ny, deltakx, deltaky
USE model, ONLY : ExBrate
IMPLICIT NONE IMPLICIT NONE
REAL(xp), INTENT(IN) :: ExBrate
! Setup the ExB shift ! Setup the ExB shearing rate and aux var
gamma_E = ExBrate
IF(gamma_E .GT. 0._xp) THEN
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)) ALLOCATE(sky_ExB(local_nky))
sky_ExB = 0._xp sky_ExB = 0._xp
! Setup the jump array and shifting flag ! Setup the jump array
ALLOCATE(jump_ExB(local_nky)) ALLOCATE(jump_ExB(local_nky))
jump_ExB = 0 jump_ExB = 0
! Setup the shifting flag array
ALLOCATE(shiftnow_ExB(local_nky)) ALLOCATE(shiftnow_ExB(local_nky))
shiftnow_ExB = .FALSE. shiftnow_ExB = .FALSE.
...@@ -39,13 +54,6 @@ CONTAINS ...@@ -39,13 +54,6 @@ CONTAINS
ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx)) ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx))
ExB_NL_factor = 1._xp ExB_NL_factor = 1._xp
inv_ExB_NL_factor = 1._xp inv_ExB_NL_factor = 1._xp
IF(ExBrate .NE. 0) THEN
t0 = deltakx/deltaky/ExBrate
inv_t0 = 1._xp/t0
ELSE ! avoid 1/0 division (t0 is killed anyway in this case)
t0 = 0._xp
inv_t0 = 0._xp
ENDIF
END SUBROUTINE Setup_ExB_shear_flow END SUBROUTINE Setup_ExB_shear_flow
...@@ -64,64 +72,61 @@ CONTAINS ...@@ -64,64 +72,61 @@ CONTAINS
IMPLICIT NONE IMPLICIT NONE
! local var ! local var
INTEGER :: iky, ikx, ikx_s, i_, loopstart, loopend, increment INTEGER :: iky, ikx, ikx_s, i_, loopstart, loopend, increment
IF(ExB) THEN
CALL start_chrono(chrono_ExBs) ! shift all fields and correct the shift value
DO iky = 1,local_Nky
! shift all fields and correct the shift value IF(shiftnow_ExB(iky)) THEN
DO iky = 1,local_Nky ! We shift the array from left to right or right to left according to the jump
IF(shiftnow_ExB(iky)) THEN ! This avoids to make copy
! We shift the array from left to right or right to left according to the jump IF(jump_ExB(iky) .GT. 0) THEN
! This avoids to make copy loopstart = 1
IF(jump_ExB(iky) .GT. 0) THEN loopend = total_nkx
loopstart = 1 increment = 1
loopend = total_nkx ELSE
increment = 1 loopstart = total_nkx
ELSE loopend = 1
loopstart = total_nkx increment = -1
loopend = 1
increment = -1
ENDIF
!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)
! so to go along the array in a monotonic way one must travel as
! 67812345 or 54321678
DO i_ = loopstart, loopend, increment
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 ENDIF
ikx_s = ikx + jump_ExB(iky) !loop to go through the array in a monotonic kx order
! We test if the shifted modes are still in contained in our resolution ! Recall: the kx array is organized as
! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN ! 6 7 8 1 2 3 4 5 (indices ikx)
IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. & ! -3 -2 -1 0 1 2 3 4 (values in dkx)
(((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. & ! so to go along the array in a monotonic way one must travel as
((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN ! 67812345 or 54321678
moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:) DO i_ = loopstart, loopend, increment
phi(iky,ikx,:) = phi(iky,ikx_s,:) IF (i_ .LT. total_nkx/2) THEN ! go to the negative kx region
psi(iky,ikx,:) = psi(iky,ikx_s,:) ikx = i_ + total_nkx/2 + 1
ELSE ! if it is not, it is lost (~dissipation for high modes) ELSE ! positive
moments(:,:,:,iky,ikx,:,:) = 0._xp ikx = i_ - total_nkx/2 + 1
phi(iky,ikx,:) = 0._xp ENDIF
psi(iky,ikx,:) = 0._xp ikx_s = ikx + jump_ExB(iky)
ENDIF ! We test if the shifted modes are still in contained in our resolution
ENDDO ! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN
! correct the shift value s(ky) for this row IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
sky_ExB(iky) = sky_ExB(iky) - jump_ExB(iky)*deltakx (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
ENDIF ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
ENDDO moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
! After shifting and correction, we update the operators and grids phi(iky,ikx,:) = phi(iky,ikx_s,:)
! update the grids psi(iky,ikx,:) = psi(iky,ikx_s,:)
CALL update_grids(sky_ExB,gxx,gxy,inv_hatB2) ELSE ! if it is not, it is lost (~dissipation for high modes)
moments(:,:,:,iky,ikx,:,:) = 0._xp
! update the EM op., the kernels and the curvature op. phi(iky,ikx,:) = 0._xp
CALL evaluate_kernels psi(iky,ikx,:) = 0._xp
CALL evaluate_EM_op ENDIF
CALL evaluate_magn_curv ENDDO
! correct the shift value s(ky) for this row
sky_ExB(iky) = sky_ExB(iky) - jump_ExB(iky)*deltakx
ENDIF
ENDDO
! After shifting and correction, we update the operators and grids
! update the grids
CALL update_grids(sky_ExB,gxx,gxy,inv_hatB2)
CALL stop_chrono(chrono_ExBs) ! update the EM op., the kernels and the curvature op.
CALL evaluate_kernels
CALL evaluate_EM_op
CALL evaluate_magn_curv
ENDIF
END SUBROUTINE Apply_ExB_shear_flow END SUBROUTINE Apply_ExB_shear_flow
! update the ExB shear value for the next time step ! update the ExB shear value for the next time step
...@@ -130,52 +135,52 @@ CONTAINS ...@@ -130,52 +135,52 @@ CONTAINS
USE grid, ONLY: local_nky, local_nky_offset, kyarray, kyarray_full, inv_dkx, xarray, Nx, Ny, & USE grid, ONLY: local_nky, local_nky_offset, kyarray, kyarray_full, inv_dkx, xarray, Nx, Ny, &
local_nx, local_nx_offset, deltax, & local_nx, local_nx_offset, deltax, &
ikyarray, inv_ikyarray, deltakx, deltaky ikyarray, inv_ikyarray, deltakx, deltaky
USE model, ONLY: ExBrate
IMPLICIT NONE IMPLICIT NONE
! local var ! local var
INTEGER :: iky, ix INTEGER :: iky, ix
REAL(xp):: dtExBshear, ky, kx, J_dp, inv_J, x REAL(xp):: dtExBshear, ky, kx, J_dp, inv_J, x
CALL start_chrono(chrono_ExBs) ! do nothing if no ExB
! update the ExB shift, jumps and flags IF(ExB) THEN
shiftnow_ExB = .FALSE. ! update the ExB shift, jumps and flags
DO iky = 1,local_Nky shiftnow_ExB = .FALSE.
sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*ExBrate*dt DO iky = 1,local_Nky
jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx) sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*gamma_E*dt
! If the jump is 1 or more for a given ky, we flag the index jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx)
! in shiftnow_ExB and will use it in Shift_fields to avoid ! If the jump is 1 or more for a given ky, we flag the index
! zero-shiftings that may be majoritary. ! in shiftnow_ExB and will use it in Shift_fields to avoid
shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0) ! zero-shiftings that may be majoritary.
ENDDO shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
ENDDO
! Update the ExB nonlinear factor... ! Update the ExB nonlinear factor...
DO iky = 1,local_Nky DO iky = 1,local_Nky
! for readability ! for readability
ky = kyarray_full(iky+local_nky_offset) ky = kyarray_full(iky+local_nky_offset)
J_dp = ikyarray(iky+local_nky_offset) J_dp = ikyarray(iky+local_nky_offset)
inv_J = inv_ikyarray(iky+local_nky_offset) inv_J = inv_ikyarray(iky+local_nky_offset)
! compute dt factor ! compute dt factor
dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp) dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
DO ix = 1,Nx DO ix = 1,Nx
x = xarray(ix) x = xarray(ix)
! assemble the ExB nonlin factor ! assemble the ExB nonlin factor
ExB_NL_factor(ix,iky) = EXP(-imagu*x*ExBrate*ky*dtExBshear) ExB_NL_factor(ix,iky) = EXP(-imagu*x*gamma_E*ky*dtExBshear)
ENDDO
ENDDO ENDDO
ENDDO ! ... and the inverse
! ... and the inverse DO iky = 1,Ny/2+1
DO iky = 1,Ny/2+1 ! for readability
! for readability ky = kyarray_full(iky)
ky = kyarray_full(iky) J_dp = ikyarray(iky)
J_dp = ikyarray(iky) inv_J = inv_ikyarray(iky)
inv_J = inv_ikyarray(iky) ! compute dt factor
! compute dt factor dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp) ky = REAL(iky-1,xp)*deltaky
ky = REAL(iky-1,xp)*deltaky DO ix = 1,local_nx
DO ix = 1,local_nx x = xarray(ix+local_nx_offset)
x = xarray(ix+local_nx_offset) ! assemble the inverse ExB nonlin factor
! assemble the inverse ExB nonlin factor inv_ExB_NL_factor(iky,ix) = EXP(imagu*x*gamma_E*ky*dtExBshear)
inv_ExB_NL_factor(iky,ix) = EXP(imagu*x*ExBrate*ky*dtExBshear) ENDDO
ENDDO ENDDO
ENDDO ENDIF
CALL stop_chrono(chrono_ExBs)
END SUBROUTINE Update_ExB_shear_flow END SUBROUTINE Update_ExB_shear_flow
END MODULE ExB_shear_flow END MODULE ExB_shear_flow
\ No newline at end of file
src/auxval.F90
+ 2
− 2
View file @ 4b80c283
...@@ -6,7 +6,7 @@ subroutine auxval ...@@ -6,7 +6,7 @@ subroutine auxval
local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset, dmax,& local_nky, local_nky_offset, total_nky, local_nkx, local_nkx_offset, dmax,&
local_nz, local_nz_offset, total_nz, init_grids_data, set_grids local_nz, local_nz_offset, total_nz, init_grids_data, set_grids
!USE array !USE array
USE model, ONLY: Na, EM, LINEARITY, N_HD USE model, ONLY: Na, EM, LINEARITY, N_HD, ExBrate
USE fourier, ONLY: init_grid_distr_and_plans USE fourier, ONLY: init_grid_distr_and_plans
use MPI, ONLY: MPI_COMM_WORLD use MPI, ONLY: MPI_COMM_WORLD
USE numerics, ONLY: build_dnjs_table, build_dv4Hp_table, compute_lin_coeff, & USE numerics, ONLY: build_dnjs_table, build_dv4Hp_table, compute_lin_coeff, &
...@@ -48,7 +48,7 @@ subroutine auxval ...@@ -48,7 +48,7 @@ subroutine auxval
! set the closure scheme in use ! set the closure scheme in use
CALL set_closure_model CALL set_closure_model
! Setup ExB shear variables ! Setup ExB shear variables
CALL Setup_ExB_shear_flow CALL Setup_ExB_shear_flow(ExBrate)
#ifdef TEST_SVD #ifdef TEST_SVD
! If we want to test SVD decomposition etc. ! If we want to test SVD decomposition etc.
CALL init_CLA(local_nky,local_np*local_nj) CALL init_CLA(local_nky,local_np*local_nj)
......
src/model_mod.F90
+ 0
− 5
View file @ 4b80c283
...@@ -31,7 +31,6 @@ MODULE model ...@@ -31,7 +31,6 @@ MODULE model
! Auxiliary variable ! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift
LOGICAL, PUBLIC, PROTECTED :: ExB = .false. ! presence of ExB background shearing rate
! Removes Landau damping in temperature and higher equation (Ivanov 2022) ! Removes Landau damping in temperature and higher equation (Ivanov 2022)
LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false. LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false.
! Flag to force the reality condition symmetry for the kx at ky=0 ! Flag to force the reality condition symmetry for the kx at ky=0
...@@ -79,10 +78,6 @@ CONTAINS ...@@ -79,10 +78,6 @@ CONTAINS
EM = .FALSE. EM = .FALSE.
ENDIF ENDIF
IF(ExBrate .GT. 0) THEN
ExB = .TRUE.
ENDIF
END SUBROUTINE model_readinputs END SUBROUTINE model_readinputs
SUBROUTINE model_outputinputs(fid) SUBROUTINE model_outputinputs(fid)
......
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