From 7387b96802b2fd42ec61502299cf42e812552891 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 7 Aug 2023 19:16:48 +0200
Subject: [PATCH] Not benchmarked ExB implementation

---
 src/ExB_shear_flow_mod.F90 | 121 ++++++++++++++++++-------------------
 1 file changed, 59 insertions(+), 62 deletions(-)

diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
index 204856da..508449fe 100644
--- a/src/ExB_shear_flow_mod.F90
+++ b/src/ExB_shear_flow_mod.F90
@@ -12,99 +12,96 @@ MODULE ExB_shear_flow
     LOGICAL,  DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift
 
     ! Routines
-    PUBLIC :: Setup_ExB_shear_flow, Update_ExB_shear_flow
+    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(g_E)
-        USE basic,      ONLY : time
-        USE grid,       ONLY : local_nky, kyarray, Lx
-        USE prec_const, ONLY : PI
+    SUBROUTINE Setup_ExB_shear_flow
+        USE grid,       ONLY : local_nky
         IMPLICIT NONE
-        REAL(xp), INTENT(IN) :: g_E ! Input shearing rate (comes from model module)
-        ! local var
-        INTEGER :: iky
-        REAL(xp):: inv_dkx
-
-        ! Store the shearing rate in this module
-        gamma_E = g_E
 
         ! Setup the ExB shift
         ALLOCATE(sky_ExB(local_nky))
         sky_ExB = 0._xp
 
         ! Setup the jump array and shifting flag
-        ALLOCATE(shiftnow_ExB(local_nky))
-        shiftnow_ExB = .FALSE.
         ALLOCATE(jump_ExB(local_nky))
         jump_ExB     = 0
+        ALLOCATE(shiftnow_ExB(local_nky))
+        shiftnow_ExB = .FALSE.
     END SUBROUTINE Setup_ExB_shear_flow
 
     ! Update according to the current ExB shear value
     ! -the grids
     ! -the spatial operators 
     ! -the fields by imposing a shift on kx
-    ! Then update the ExB shear value for the next time step
-    SUBROUTINE Update_ExB_shear_flow
-        USE basic,      ONLY : dt
-        USE grid,       ONLY : local_nky, kyarray, update_grids
-        USE prec_const, ONLY : PI
-        USE geometry,   ONLY : gxx,gxy,gyy,inv_hatB2, evaluate_magn_curv
-        USE numerics,   ONLY : evaluate_EM_op, evaluate_kernels
+    SUBROUTINE Apply_ExB_shear_flow
+        USE basic,      ONLY: chrono_ExBs, start_chrono, stop_chrono
+        USE grid,       ONLY: local_nky, kxarray, update_grids, &
+            local_nkx, deltakx, kx_min, kx_max
+        USE prec_const, ONLY: PI
+        USE geometry,   ONLY: gxx,gxy,inv_hatB2, evaluate_magn_curv
+        USE numerics,   ONLY: evaluate_EM_op, evaluate_kernels
+        USE fields,     ONLY: moments, phi, psi
         IMPLICIT NONE
         ! local var
-        INTEGER :: iky
-        REAL(xp):: inv_dkx
+        INTEGER :: iky, ikx, ikx_s
 
-        ! update the grids
-        CALL update_grids(sky_ExB(iky),gxx,gxy,gyy,inv_hatB2)
+        CALL start_chrono(chrono_ExBs)
 
-        ! update the EM operators and the kernels
+        ! shift all fields and correct the shift value
+        DO iky = 1,local_Nky
+            IF(shiftnow_ExB(iky)) THEN
+                print*, "SHIFT ARRAYS"
+                ! shift all fields
+                DO ikx = 1,local_nkx
+                    ikx_s = ikx + jump_ExB(iky)
+                    IF( (kxarray(iky,ikx) .GE. kx_min) .AND. (kxarray(iky,ikx) .LE. kx_max) ) THEN
+                        moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
+                        phi(iky,ikx,:)             = phi(iky,ikx_s,:)
+                        psi(iky,ikx,:)             = psi(iky,ikx_s,:)
+                    ELSE
+                        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
 
-        ! shift all fields and correct the shift value
-        CALL Shift_fields
+        CALL stop_chrono(chrono_ExBs)
+    END SUBROUTINE Apply_ExB_shear_flow
 
-        ! update the ExB shift
+    ! update the ExB shear value for the next time step
+    SUBROUTINE Update_ExB_shear_flow
+        USE basic,      ONLY: dt, chrono_ExBs, start_chrono, stop_chrono
+        USE grid,       ONLY: local_nky, kyarray, inv_dkx
+        USE model,      ONLY: ExBrate
+        IMPLICIT NONE
+        ! local var
+        INTEGER :: iky
+        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)*gamma_E*dt
+            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
-        CONTAINS
-
-        SUBROUTINE Shift_fields
-            USE grid,  ONLY: local_nky, kxarray, kyarray, local_nkx, deltakx, &
-                             kx_min, kx_max
-            USE fields,ONLY: moments, phi, psi
-            IMPLICIT NONE
-            ! local var
-            INTEGER :: iky, ikx, ikx_s
-            REAL(xp):: inv_dkx
-            DO iky = 1,local_Nky
-                IF(shiftnow_ExB(iky)) THEN
-                    ! shift all fields
-                    DO ikx = 1,local_nkx
-                        ikx_s = ikx + jump_ExB(iky)
-                        IF( (kxarray(iky,ikx) .GE. kx_min) .AND. (kxarray(iky,ikx) .LE. kx_max) ) THEN
-                            moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
-                            phi(iky,ikx,:)             = phi(iky,ikx_s,:)
-                            psi(iky,ikx,:)             = psi(iky,ikx_s,:)
-                        ELSE
-                            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
-        END SUBROUTINE Shift_fields
+        CALL stop_chrono(chrono_ExBs)
     END SUBROUTINE Update_ExB_shear_flow
-
-
 END MODULE ExB_shear_flow
\ No newline at end of file
-- 
GitLab