From 4b80c28354e9828f457d656c75f9876aa4b1d31f Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 13 Sep 2023 09:59:38 +0200
Subject: [PATCH] move the ExB flag
---
src/ExB_shear_flow_mod.F90 | 221 +++++++++++++++++++------------------
src/auxval.F90 | 4 +-
src/model_mod.F90 | 5 -
3 files changed, 115 insertions(+), 115 deletions(-)
diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 7cc0033a..5236d28e 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -13,24 +13,39 @@ MODULE ExB_shear_flow
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
! Routines
PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
CONTAINS
! 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 model, ONLY : ExBrate
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))
sky_ExB = 0._xp
- ! Setup the jump array and shifting flag
+ ! Setup the jump array
ALLOCATE(jump_ExB(local_nky))
- jump_ExB = 0
+ jump_ExB = 0
+
+ ! Setup the shifting flag array
ALLOCATE(shiftnow_ExB(local_nky))
shiftnow_ExB = .FALSE.
@@ -39,13 +54,6 @@ CONTAINS
ALLOCATE(inv_ExB_NL_factor(Ny/2+1,local_nx))
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
@@ -64,64 +72,61 @@ CONTAINS
IMPLICIT NONE
! local var
INTEGER :: iky, ikx, ikx_s, i_, loopstart, loopend, increment
-
- CALL start_chrono(chrono_ExBs)
-
- ! shift all fields and correct the 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
- ! This avoids to make copy
- IF(jump_ExB(iky) .GT. 0) THEN
- loopstart = 1
- loopend = total_nkx
- increment = 1
- ELSE
- loopstart = total_nkx
- 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
+ IF(ExB) THEN
+ ! shift all fields and correct the 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
+ ! This avoids to make copy
+ IF(jump_ExB(iky) .GT. 0) THEN
+ loopstart = 1
+ loopend = total_nkx
+ increment = 1
+ ELSE
+ loopstart = total_nkx
+ loopend = 1
+ increment = -1
ENDIF
- ikx_s = ikx + jump_ExB(iky)
- ! We test if the shifted modes are still in contained in our resolution
- ! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN
- 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
- moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
- phi(iky,ikx,:) = phi(iky,ikx_s,:)
- psi(iky,ikx,:) = psi(iky,ikx_s,:)
- ELSE ! if it is not, it is lost (~dissipation for high modes)
- moments(:,:,:,iky,ikx,:,:) = 0._xp
- phi(iky,ikx,:) = 0._xp
- psi(iky,ikx,:) = 0._xp
- ENDIF
- 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)
-
- ! update the EM op., the kernels and the curvature op.
- CALL evaluate_kernels
- CALL evaluate_EM_op
- CALL evaluate_magn_curv
+ !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
+ ikx_s = ikx + jump_ExB(iky)
+ ! We test if the shifted modes are still in contained in our resolution
+ ! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN
+ 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
+ moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
+ phi(iky,ikx,:) = phi(iky,ikx_s,:)
+ psi(iky,ikx,:) = psi(iky,ikx_s,:)
+ ELSE ! if it is not, it is lost (~dissipation for high modes)
+ moments(:,:,:,iky,ikx,:,:) = 0._xp
+ phi(iky,ikx,:) = 0._xp
+ psi(iky,ikx,:) = 0._xp
+ ENDIF
+ 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
! update the ExB shear value for the next time step
@@ -130,52 +135,52 @@ CONTAINS
USE grid, ONLY: local_nky, local_nky_offset, kyarray, kyarray_full, inv_dkx, xarray, Nx, Ny, &
local_nx, local_nx_offset, deltax, &
ikyarray, inv_ikyarray, deltakx, deltaky
- USE model, ONLY: ExBrate
IMPLICIT NONE
! local var
INTEGER :: iky, ix
REAL(xp):: dtExBshear, ky, kx, J_dp, inv_J, x
- CALL start_chrono(chrono_ExBs)
- ! update the ExB shift, jumps and flags
- shiftnow_ExB = .FALSE.
- DO iky = 1,local_Nky
- sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*ExBrate*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
+ ! do nothing if no ExB
+ IF(ExB) THEN
+ ! update 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)
+ ENDDO
- ! Update the ExB nonlinear factor...
- DO iky = 1,local_Nky
- ! for readability
- ky = kyarray_full(iky+local_nky_offset)
- J_dp = ikyarray(iky+local_nky_offset)
- inv_J = inv_ikyarray(iky+local_nky_offset)
- ! compute dt factor
- dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
- DO ix = 1,Nx
- x = xarray(ix)
- ! assemble the ExB nonlin factor
- ExB_NL_factor(ix,iky) = EXP(-imagu*x*ExBrate*ky*dtExBshear)
+ ! Update the ExB nonlinear factor...
+ DO iky = 1,local_Nky
+ ! for readability
+ ky = kyarray_full(iky+local_nky_offset)
+ J_dp = ikyarray(iky+local_nky_offset)
+ inv_J = inv_ikyarray(iky+local_nky_offset)
+ ! compute dt factor
+ dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
+ DO ix = 1,Nx
+ x = xarray(ix)
+ ! assemble the ExB nonlin factor
+ ExB_NL_factor(ix,iky) = EXP(-imagu*x*gamma_E*ky*dtExBshear)
+ ENDDO
ENDDO
- ENDDO
- ! ... and the inverse
- DO iky = 1,Ny/2+1
- ! for readability
- ky = kyarray_full(iky)
- J_dp = ikyarray(iky)
- inv_J = inv_ikyarray(iky)
- ! compute dt factor
- dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
- ky = REAL(iky-1,xp)*deltaky
- 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*ExBrate*ky*dtExBshear)
+ ! ... and the inverse
+ DO iky = 1,Ny/2+1
+ ! for readability
+ ky = kyarray_full(iky)
+ J_dp = ikyarray(iky)
+ inv_J = inv_ikyarray(iky)
+ ! compute dt factor
+ dtExBshear = time - t0*inv_J*ANINT(J_dp*time*inv_t0,xp)
+ ky = REAL(iky-1,xp)*deltaky
+ 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*ky*dtExBshear)
+ ENDDO
ENDDO
- ENDDO
- CALL stop_chrono(chrono_ExBs)
+ ENDIF
END SUBROUTINE Update_ExB_shear_flow
END MODULE ExB_shear_flow
\ No newline at end of file
diff --git a/src/auxval.F90 b/src/auxval.F90
index 32966633..0cf3960c 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -6,7 +6,7 @@ subroutine auxval
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
!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 MPI, ONLY: MPI_COMM_WORLD
USE numerics, ONLY: build_dnjs_table, build_dv4Hp_table, compute_lin_coeff, &
@@ -48,7 +48,7 @@ subroutine auxval
! set the closure scheme in use
CALL set_closure_model
! Setup ExB shear variables
- CALL Setup_ExB_shear_flow
+ CALL Setup_ExB_shear_flow(ExBrate)
#ifdef TEST_SVD
! If we want to test SVD decomposition etc.
CALL init_CLA(local_nky,local_np*local_nj)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 4072421a..f468a0cd 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -31,7 +31,6 @@ MODULE model
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
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)
LOGICAL, PUBLIC, PROTECTED :: RM_LD_T_EQ = .false.
! Flag to force the reality condition symmetry for the kx at ky=0
@@ -79,10 +78,6 @@ CONTAINS
EM = .FALSE.
ENDIF
- IF(ExBrate .GT. 0) THEN
- ExB = .TRUE.
- ENDIF
-
END SUBROUTINE model_readinputs
SUBROUTINE model_outputinputs(fid)
--
GitLab