epochpic · usinhampik · May 8, 2020 · May 8, 2020 · May 26, 2020 · May 28, 2020
diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
@@ -14,3 +14,5 @@ Phil Tooley <phil.tooley@gmail.com>
 Stephan Kuschel <stephan.kuschel@gmail.com>
 Stuart Morris <sjm630@york.ac.uk>
 Tom Goffrey <t.goffrey@warwick.ac.uk>
+Ujjwal Sinha <u.sinha@fz-juelich.de> (Acknowledgement:
+Contributions from Ujjwal Sinha was financially supported by the PRACE project funded in part by the EU’s Horizon 2020 Research and Innovation programme (2014-2020) under grant agreement 823767.)
diff --git a/README.md b/README.md
@@ -1,3 +1,19 @@
+## EPOCH development
+
+**This is not the official repository for the EPOCH code.**
+
+This version of EPOCH is being developed within the project PiCKeX-
+Particle Kinetic codes for Exascale plasma simulation, which is a
+part of Work Package 8 (WP8) of PRACE-6IP. This repository has been
+cloned to the Jülich Supercomputing Centre (JSC) GitLab server to
+perform and publish this work. The main production version of EPOCH is
+maintained by the University of Warwick. The development version found
+here might not have the full functionality of the production version of
+EPOCH, but includes a performance improved, refactored branch that will
+be merged with the official version by the end of the project.
+
+Current developments include improving the efficiency of the moving window book-keeping, reducing the number MPI calls for the field boundary conditions and the particle boundary conditions.
+
 # *** PLEASE READ THIS NOTE ***
 
 If you are obtaining this code from the github server *DO NOT* use the

diff --git a/epoch2d/Makefile b/epoch2d/Makefile
@@ -50,11 +50,11 @@ endif
 # Intel
 # =====
 ifeq ($(strip $(COMPILER)),intel)
-  FFLAGS = -O3 -g -stand f03
+  FFLAGS = -O3 -g -stand f03 -qopenmp
   ifeq ($(strip $(CONS)),1)
     FLTCONS = -pc64 -fltconsistency
   endif
-  #FFLAGS = -O3 -heap-arrays 64 -ipo -xHost # Optimised (B)
+  #FFLAGS = -O3 -heap-arrays 64  -qopt-report-phase=vec -qopt-report=5 -qopenmp # Optimised (B)
   #FFLAGS = -O3 -heap-arrays 64 -ipo -xAVX  # Optimised (W)
   ifeq ($(strip $(MODE)),debug)
     FFLAGS = -O0 -fpe0 -nothreads -traceback -fltconsistency \
@@ -554,4 +554,4 @@ tokenizer_blocks.o: tokenizer_blocks.f90 strings.o
 utilities.o: utilities.f90 constants.o shared_data.o
 version_data.o: version_data.F90 $(SRCDIR)/COMMIT
 welcome.o: welcome.F90 shared_data.o terminal_controls.o version_data.o
-window.o: window.F90 boundary.o evaluate.o partlist.o
+window.o: window.F90 boundary.o evaluate.o partlist.o deck_io_block.o
diff --git a/epoch2d/src/boundary.F90 b/epoch2d/src/boundary.F90
@@ -316,6 +316,218 @@ END SUBROUTINE do_field_mpi_with_lengths
 
 
 
+  SUBROUTINE all_comp_field_bc(fieldx, fieldy, fieldz, ng, nx_local, &
+        ny_local)
+
+
+    !            |     field_top      |           
+    !____________|____________________|____________
+    !            |                    |
+    ! field_left |                    | field_right
+    !____________|____________________|____________
+    !            |                    |
+    !            |     field_bottom   |   
+
+    INTEGER, INTENT(IN) :: ng
+    REAL(num), DIMENSION(1-ng:,1-ng:), INTENT(INOUT) :: fieldx, fieldy, fieldz
+    INTEGER, INTENT(IN) :: nx_local, ny_local
+    INTEGER, DIMENSION(c_ndims) :: sizes, subsizes
+    INTEGER :: basetype, sz, szmax, i, j, k, n
+    REAL(num), ALLOCATABLE :: field(:)
+    REAL(num), ALLOCATABLE :: temp(:)
+    INTEGER :: xmin, xmax, ymin, ymax, offset
+
+    basetype = mpireal
+
+    sizes(1) = nx_local + 2 * ng
+    sizes(2) = ny_local + 2 * ng
+
+    szmax = 3 * sizes(1) * ng
+    sz = 3 * sizes(2) * ng
+    IF (sz > szmax) szmax = sz
+
+    ALLOCATE(temp(szmax))
+    ALLOCATE(field(szmax))
+
+    subsizes(1) = ng
+    subsizes(2) = sizes(2)
+
+    sz = 3 * subsizes(1) * subsizes(2)
+
+!    offset0 = 0
+!    offset1 = subsizes(1) * subsizes(2)
+!    offset2 = 2 * offset1
+
+    offset = subsizes(1) * subsizes(2)
+
+    xmin = 1
+    xmax = ng
+    ymin = 1-ng
+    ymax = subsizes(2)-ng
+
+    CALL load_field_boundaries_to_buffer(fieldx, field, &
+        xmin, xmax, ymin, ymax, 0*offset)
+    CALL load_field_boundaries_to_buffer(fieldy, field, &
+        xmin, xmax, ymin, ymax, offset)
+    CALL load_field_boundaries_to_buffer(fieldz, field, &
+        xmin, xmax, ymin, ymax, 2*offset)
+
+    CALL MPI_SENDRECV(field, sz, basetype, proc_x_min, &
+        tag, temp, sz, basetype, proc_x_max, tag, comm, status, errcode)
+
+    xmin = nx_local + 1
+    xmax = subsizes(1) + nx_local
+
+    IF (.NOT. x_max_boundary .OR. bc_field(c_bd_x_max)==c_bc_periodic) THEN
+
+      CALL unload_field_boundaries_from_buffer(fieldx, temp, &
+          xmin, xmax, ymin, ymax, 0*offset)
+      CALL unload_field_boundaries_from_buffer(fieldy, temp, &
+          xmin, xmax, ymin, ymax, offset)
+      CALL unload_field_boundaries_from_buffer(fieldz, temp, &
+          xmin, xmax, ymin, ymax, 2*offset)
+
+    END IF
+
+    xmin = nx_local - ng + 1
+    xmax = nx_local
+
+    CALL load_field_boundaries_to_buffer(fieldx, field, &
+        xmin, xmax, ymin, ymax, 0*offset)
+    CALL load_field_boundaries_to_buffer(fieldy, field, &
+        xmin, xmax, ymin, ymax, offset)
+    CALL load_field_boundaries_to_buffer(fieldz, field, &
+        xmin, xmax, ymin, ymax, 2*offset)
+
+    CALL MPI_SENDRECV(field, sz, basetype, proc_x_max, &
+        tag, temp, sz, basetype, proc_x_min, tag, comm, status, errcode)
+
+    xmin = 1-ng
+    xmax = subsizes(1)-ng
+
+    IF (.NOT. x_min_boundary .OR. bc_field(c_bd_x_min)==c_bc_periodic) THEN
+
+      CALL unload_field_boundaries_from_buffer(fieldx, temp, &
+          xmin, xmax, ymin, ymax, 0*offset)
+      CALL unload_field_boundaries_from_buffer(fieldy, temp, &
+          xmin, xmax, ymin, ymax, offset)
+      CALL unload_field_boundaries_from_buffer(fieldz, temp, &
+        xmin, xmax, ymin, ymax, 2*offset)
+
+    END IF
+
+    subsizes(1) = sizes(1)
+    subsizes(2) = ng
+    sz = 3 * subsizes(1) * subsizes(2)
+
+!    offset0 = 0
+!    offset1 = subsizes(1) * subsizes(2)
+!    offset2 = 2 * offset1
+
+    offset = subsizes(1) * subsizes(2)   
+
+    xmin = 1 - ng
+    xmax = subsizes(1) - ng 
+    ymin = ny_local - subsizes(2) + 1
+    ymax = ny_local
+
+    CALL load_field_boundaries_to_buffer(fieldx, field, &
+        xmin, xmax, ymin, ymax, 0*offset)
+    CALL load_field_boundaries_to_buffer(fieldy, field, &
+        xmin, xmax, ymin, ymax, offset)
+    CALL load_field_boundaries_to_buffer(fieldz, field, &
+        xmin, xmax, ymin, ymax, 2*offset)
+
+    CALL MPI_SENDRECV(field, sz, basetype, proc_y_max, &
+        tag, temp, sz, basetype, proc_y_min, tag, comm, status, errcode)
+
+    ymin = 1 - ng
+    ymax = subsizes(2) - ng
+
+    IF (.NOT. y_min_boundary .OR. bc_field(c_bd_y_min)==c_bc_periodic) THEN
+
+      CALL unload_field_boundaries_from_buffer(fieldx, temp, &
+          xmin, xmax, ymin, ymax, 0*offset)
+      CALL unload_field_boundaries_from_buffer(fieldy, temp, &
+          xmin, xmax, ymin, ymax, offset)
+      CALL unload_field_boundaries_from_buffer(fieldz, temp, &
+          xmin, xmax, ymin, ymax, 2*offset)
+
+    END IF
+
+    ymin = 1
+    ymax = subsizes(2)
+
+    CALL load_field_boundaries_to_buffer(fieldx, field, &
+        xmin, xmax, ymin, ymax, 0*offset)
+    CALL load_field_boundaries_to_buffer(fieldy, field, &
+        xmin, xmax, ymin, ymax, offset)
+    CALL load_field_boundaries_to_buffer(fieldz, field, &
+        xmin, xmax, ymin, ymax, 2*offset)
+
+    ymin = ny_local + 1
+    ymax = subsizes(2) + ny_local
+
+    CALL MPI_SENDRECV(field, sz, basetype, proc_y_min, &
+        tag, temp, sz, basetype, proc_y_max, tag, comm, status, errcode)
+
+    IF (.NOT. y_max_boundary .OR. bc_field(c_bd_y_max)==c_bc_periodic) THEN
+
+      CALL unload_field_boundaries_from_buffer(fieldx, temp, &
+          xmin, xmax, ymin, ymax, 0*offset)
+      CALL unload_field_boundaries_from_buffer(fieldy, temp, &
+          xmin, xmax, ymin, ymax, offset)
+      CALL unload_field_boundaries_from_buffer(fieldz, temp, &
+          xmin, xmax, ymin, ymax, 2*offset)
+
+    END IF
+
+
+    DEALLOCATE(field)
+    DEALLOCATE(temp)
+
+  END SUBROUTINE all_comp_field_bc
+
+
+
+  SUBROUTINE load_field_boundaries_to_buffer(field, buffer, &
+            xmin, xmax, ymin, ymax, offset)
+
+    INTEGER, INTENT(IN) :: xmin, xmax, ymin, ymax, offset
+    REAL(num), DIMENSION(1-ng:, 1-ng:), INTENT(INOUT) :: field
+    REAL(num), DIMENSION(:), INTENT(INOUT) :: buffer
+    INTEGER :: i, j, n
+
+    n = 1
+    DO j = ymin, ymax
+    DO i = xmin, xmax
+      buffer(n + offset) = field(i,j)
+      n = n + 1
+    END DO
+    END DO
+
+   END SUBROUTINE load_field_boundaries_to_buffer
+
+   SUBROUTINE unload_field_boundaries_from_buffer(field, buffer, &
+            xmin, xmax, ymin, ymax, offset)
+
+     INTEGER, INTENT(IN) :: xmin, xmax, ymin, ymax, offset
+     REAL(num), DIMENSION(1-ng:, 1-ng:), INTENT(INOUT) :: field
+     REAL(num), DIMENSION(:), INTENT(INOUT) :: buffer
+     INTEGER :: i, j, n
+
+     n = 1
+     DO j = ymin, ymax
+     DO i = xmin, xmax
+       field(i,j) = buffer(n + offset)
+       n = n + 1
+     END DO
+     END DO
+
+   END SUBROUTINE unload_field_boundaries_from_buffer
+
+
+
   SUBROUTINE do_field_mpi_with_lengths_r4(field, ng, nx_local, ny_local)
 
     INTEGER, INTENT(IN) :: ng
@@ -810,9 +1022,10 @@ SUBROUTINE efield_bcs
     INTEGER :: i
 
     ! These are the MPI boundaries
-    CALL field_bc(ex, ng)
-    CALL field_bc(ey, ng)
-    CALL field_bc(ez, ng)
+    !CALL field_bc(ex, ng)
+    !CALL field_bc(ey, ng)
+    !CALL field_bc(ez, ng)
+    CALL all_comp_field_bc(ex, ey, ez, ng, nx, ny)
 
     ! Perfectly conducting boundaries
     DO i = c_bd_x_min, c_bd_x_max, c_bd_x_max - c_bd_x_min
@@ -861,9 +1074,11 @@ SUBROUTINE bfield_bcs(mpi_only)
     INTEGER :: i
 
     ! These are the MPI boundaries
-    CALL field_bc(bx, ng)
-    CALL field_bc(by, ng)
-    CALL field_bc(bz, ng)
+    !CALL field_bc(bx, ng)
+    !CALL field_bc(by, ng)
+    !CALL field_bc(bz, ng)
+    CALL all_comp_field_bc(bx, by, bz, ng, nx, ny)
+
 
     IF (mpi_only) RETURN
 

diff --git a/epoch2d/src/deck/deck_io_block.F90 b/epoch2d/src/deck/deck_io_block.F90
@@ -1010,13 +1010,15 @@ SUBROUTINE init_io_block(io_block)
     io_block%name = ''
     io_block%dt_snapshot = -1.0_num
     io_block%time_prev = 0.0_num
+    io_block%buffer_time_prev = 0.0_num
     io_block%time_first = 0.0_num
     io_block%dt_average = -1.0_num
     io_block%dt_min_average = -1.0_num
     io_block%average_time = -1.0_num
     io_block%average_time_start = -1.0_num
     io_block%nstep_snapshot = -1
     io_block%nstep_prev = 0
+    io_block%buffer_nstep_prev = 0
     io_block%nstep_first = 0
     io_block%nstep_average = -1
     io_block%restart = .FALSE.

diff --git a/epoch2d/src/deck/strings.f90 b/epoch2d/src/deck/strings.f90
@@ -323,17 +323,15 @@ END FUNCTION lowercase
 
 
 
-  FUNCTION trim_string(string)
+  FUNCTION trim_string(string_in)
 
     CHARACTER(LEN=c_max_string_length) :: trim_string
-    CHARACTER(LEN=*) :: string
-    CHARACTER(LEN=:), ALLOCATABLE :: string_copy
+    CHARACTER(LEN=*) :: string_in
+    CHARACTER(LEN=LEN_TRIM(string_in)) :: string
 
-    ALLOCATE(CHARACTER(LEN=LEN(string)) :: string_copy)
-
-    string_copy = ADJUSTL(string)
-    IF (LEN_TRIM(string_copy) > c_max_string_length) THEN
-      trim_string = string_copy(1:c_max_string_length)
+    string = ADJUSTL(string_in)
+    IF (LEN_TRIM(string) > c_max_string_length) THEN
+      trim_string = string(1:c_max_string_length)
     ELSE
       trim_string = TRIM(string_copy)
     END IF

diff --git a/epoch2d/src/epoch2d.F90 b/epoch2d/src/epoch2d.F90
@@ -152,7 +152,7 @@ PROGRAM pic
     ELSE
       time = time + dt / 2.0_num
       CALL update_eb_fields_final
-      CALL moving_window
+      CALL moving_window(step)
     END IF
   ELSE
     dt_store = dt
@@ -258,7 +258,8 @@ PROGRAM pic
 
     CALL update_eb_fields_final
 
-    CALL moving_window
+    CALL moving_window(step)
+
   END DO
 
   IF (rank == 0) runtime = MPI_WTIME() - walltime_started

diff --git a/epoch2d/src/housekeeping/current_smooth.F90 b/epoch2d/src/housekeeping/current_smooth.F90
@@ -30,9 +30,10 @@ SUBROUTINE current_finish
 
     CALL current_bcs
 
-    CALL field_bc(jx, jng)
-    CALL field_bc(jy, jng)
-    CALL field_bc(jz, jng)
+    !CALL field_bc(jx, jng)
+    !CALL field_bc(jy, jng)
+    !CALL field_bc(jz, jng)
+    CALL all_comp_field_bc(jx, jy, jz, jng, nx, ny)
 
     IF (smooth_currents) CALL smooth_current
 

diff --git a/epoch2d/src/housekeeping/setup.F90 b/epoch2d/src/housekeeping/setup.F90
@@ -990,13 +990,17 @@ SUBROUTINE restart_data(step)
       END IF
       IF (io_block_list(i)%dt_snapshot > 0.0_num) THEN
         io_block_list(i)%time_prev = time
+        io_block_list(i)%buffer_time_prev = time
       ELSE
         io_block_list(i)%time_prev = 0.0_num
+        io_block_list(i)%buffer_time_prev = 0.0_num
       END IF
       IF (io_block_list(i)%nstep_snapshot > 0) THEN
         io_block_list(i)%nstep_prev = step
+        io_block_list(i)%buffer_nstep_prev = step
       ELSE
         io_block_list(i)%nstep_prev = 0
+        io_block_list(i)%buffer_nstep_prev = 0
       END IF
       io_block_list(i)%walltime_prev = time
       IF (ALLOCATED(io_block_list(i)%dump_at_nsteps)) THEN