Changeset 1228ef7 in flexpart.git

Timestamp:

Aug 2, 2021, 12:11:23 PM (3 years ago)

Author:

Espen Sollum <eso@…>

Branches:

dev

Parents:

6cb0801

Message:

MPI: fix for mquasilag output

Location:

src

Files:

• makefile (modified) (2 diffs)
• mpi_mod.f90 (modified) (4 diffs)
• partoutput_short_mpi.f90 (modified) (4 diffs)
• releaseparticles_mpi.f90 (modified) (4 diffs)
• timemanager_mpi.f90 (modified) (2 diffs)

src/makefile

@@ -90 +90 @@
 ## OPTIMIZATION LEVEL
 O_LEV = 2 # [0,1,2,3,g,s,fast]
-O_LEV_DBG = g # [0,g]
+O_LEV_DBG = 0 # [0,g]
 
 ## LIBRARIES
@@ -100 +100 @@
 
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
-LDDEBUG  = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
+LDDEBUG  = $(DBGFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
 
 MODOBJS = \

src/mpi_mod.f90

@@ -178 +178 @@
     implicit none
 
-    integer :: i,j,s,addmaxpart=0
+    integer :: i,j,s
 
 ! Each process gets an ID (mp_pid) in the range 0,..,mp_np-1
@@ -486 +486 @@
 
     integer,intent(in) :: num_part
-    integer :: i,jj, addone
+    integer :: i
 
 ! Exit if too many particles
@@ -589 +589 @@
     integer, intent(in) :: itime
     real :: pmin,z
-    integer :: i,jj,nn, num_part=1,m,imin, num_trans
-    logical :: first_iter
+    integer :: i,jj,nn,m,imin, num_trans
     integer,allocatable,dimension(:) :: idx_arr
     real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
@@ -613 +612 @@
     end do
 
-! Do not rebalance particles for ipout=3    
-    if (ipout.eq.3) return
+! Do not rebalance particles for ipout=3 or mquasilag=1
+    if ((ipout.eq.3).or.(mquasilag.eq.1)) return
 
 ! For each successive element in index array, see if a lower value exists

src/partoutput_short_mpi.f90

@@ -13 +13 @@
   !                                                                            *
   !     12/2014 eso: Version for MPI                                           *
-  !                  Processes sequentially access and append data to file     *
-  !                  NB: Do not use yet!                                       *
+  !                  Particle positions are sent to root process for output    *
+  !                                                                            *
   !*****************************************************************************
   !                                                                            *
@@ -28 +28 @@
 
   real(kind=dp) :: jul
-  integer :: itime,i,j,jjjjmmdd,ihmmss,numshortout,numshortall
+  integer, dimension(:), allocatable :: numshorts,displs
+  integer :: itime,i,j,jjjjmmdd,ihmmss,numshortout,numshortall,numshortmpi
   integer :: ix,jy,ixp,jyp
   real :: xlon,ylat,zlim,dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,topo
@@ -35 +36 @@
   integer(kind=2) :: idump(3,maxpart)
   integer :: i4dump(maxpart)
-  character(LEN=8) :: file_stat='OLD'
+  integer(kind=2),dimension(:,:),allocatable :: idump_all(:,:)
+  integer,dimension(:), allocatable :: i4dump_all(:)
+!  character(LEN=8) :: file_stat='OLD'
+  character(LEN=8) :: file_stat='REPLACE'
 
-  ! MPI root process creates the file, other processes append to it
-  if (lroot) file_stat='REPLACE'
+! This is not needed, in this version only root process writes the file
+  ! if (lroot) then
+  !   file_stat='REPLACE'
+  ! end if
+
+! Array to gather numshortout from all processes
+  allocate(numshorts(mp_partgroup_np), displs(mp_partgroup_np))
 
   ! Determine current calendar date, needed for the file name
@@ -115 +124 @@
         i4dump(numshortout)=npoint(i)
       endif
-
     endif
   end do
 
 
+  
+! Get total number of particles from all processes
+!************************************************
+  call MPI_Allgather(numshortout, 1, MPI_INTEGER, numshorts, 1, MPI_INTEGER, &
+       mp_comm_used, mp_ierr)
+  
+  numshortmpi = sum(numshorts(:))
+
+
+! Gather all data at root process
+!********************************
+  allocate(idump_all(3,numshortmpi), i4dump_all(numshortmpi))
+  displs(1)=0
+  do i=2,mp_partgroup_np
+    displs(i)=displs(i-1)+numshorts(i-1)
+  end do
+
+  call MPI_gatherv(i4dump, numshortout, MPI_INTEGER, i4dump_all, numshorts(:), &
+       & displs, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
+  displs = displs*3  
+  call MPI_gatherv(idump, 3*numshortout, MPI_INTEGER2, idump_all, 3*numshorts(:), &
+       & displs, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
+  
   ! Open output file and write the output
   !**************************************
 
-  open(unitshortpart,file=path(2)(1:length(2))//'shortposit_'//adate// &
-       atime,form='unformatted',status=file_stat,position='append')
+  if (lroot) then ! MPI root process only
+    open(unitshortpart,file=path(2)(1:length(2))//'shortposit_'//adate// &
+         atime,form='unformatted',status=file_stat,position='append')
+    write(unitshortpart) itime 
+    write(unitshortpart) numshortmpi
+    write(unitshortpart) &
+           (i4dump_all(i),(idump_all(j,i),j=1,3),i=1,numshortmpi)
+    close(unitshortpart)
+  end if
 
-  ! Write current time to file
-  !***************************
-
-  if (lroot) write(unitshortpart) itime ! MPI root process only
-  ! :TODO: get total numshortout (MPI reduction), add MPI barrier, open file
-  ! sequentially below
-  write(unitshortpart) numshortout 
-  write(unitshortpart) &
-       (i4dump(i),(idump(j,i),j=1,3),i=1,numshortout)
-
-
-  write(*,*) numshortout,numshortall
-
-  close(unitshortpart)
+  deallocate(idump_all, i4dump_all)
 
 end subroutine partoutput_short

src/releaseparticles_mpi.f90

@@ -51 +51 @@
   real(kind=dp) :: juldate,julmonday,jul,jullocal,juldiff
   real,parameter :: eps=nxmax/3.e5,eps2=1.e-6
-  integer :: mind2
+  integer :: mind2, numpartcount_mpi
 ! mind2        eso: pointer to 2nd windfield in memory
-
   integer :: idummy = -7
 !save idummy,xmasssave
@@ -60 +59 @@
   logical :: first_call=.true.
 
-! Use different seed for each process.
-!****************************************************************************
+! Use different random seed for each process
+!*******************************************
   if (first_call) then
     idummy=idummy+mp_seed
@@ -68 +67 @@
 
   mind2=memind(2)
+
+! For mquasilag=1, assign unique particle ID across processes
+  numpartcount_mpi=mp_partid-mp_partgroup_np+1
 
 ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time)
@@ -210 +212 @@
                  nclassunc)
             numparticlecount=numparticlecount+1
+! Use a stride equal to number of processes for the MPI version
+            numpartcount_mpi=numpartcount_mpi+mp_partgroup_np
             if (mquasilag.eq.0) then
               npoint(ipart)=i
             else
-              npoint(ipart)=numparticlecount
+              npoint(ipart)=numpartcount_mpi
             endif
             idt(ipart)=mintime               ! first time step

src/timemanager_mpi.f90

@@ -448 +448 @@
       endif
 
-! :TODO: MPI output of particle positions;  each process sequentially
-!   access the same file
       if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
-           call partoutput_short(itime)    ! dump particle positions in extremely compressed format
+           call partoutput_short(itime) ! dump particle positions in extremely compressed format
 
 
@@ -557 +555 @@
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
         if ((iflux.eq.1).and.(lnetcdfout.eq.0)) call fluxoutput(itime)
+#ifdef USE_NCF
         if ((iflux.eq.1).and.(lnetcdfout.eq.1)) call fluxoutput_netcdf(itime)
+#endif
         if (mp_measure_time) call mpif_mtime('iotime',1)