!================================================================================ module m_particles !================================================================================ ! Incompressible pseudo-spectral code with scalars ! ! The module contains subroutines for Largangian particle tracking. ! ! !================================================================================ ! VARIABLES !================================================================================ use m_parameters use m_io implicit none ! logical variable that indicates if we integrate the particles ! logical :: int_particles ! when to release the particles real*8 :: time_p ! array which contains the particles' positions in the grid units real*8, allocatable :: xyzp(:,:) ! array with second-layer velocities for particles real*8, allocatable :: uvwp(:,:) ! array with particles' tags (1 to nptot) integer*4, allocatable :: ipart(:) ! array with particles' addresses (which processors hold the particle) integer*4, allocatable :: myid_part(:), itmp_part(:) ! min and max z-coordinate for a slab real*8 :: zmin,zmax ! temp arrays for particle interpolation real*8, allocatable :: wrk1p(:,:), wrk2p(:,:), wrk3p(:,:) ! tmp arrays for CINT interpolation procedure (see particles_move_cint.f90) ! real*8, allocatable :: ! needed to avoid deadlock in send/receive particles real*8, allocatable :: xyzp1(:,:), uvwp1(:,:) integer*4, allocatable :: ipart1(:) ! number of particles to be sent/received to/from neighbours when moving ! the particles integer*4 :: np_send_u, np_send_d, np_get_u, np_get_d ! the last iteration at which the particles were written out integer :: particles_last_dump = -1 !================================================================================ !================================================================================ contains !================================================================================ ! SUBROUTINES !================================================================================ subroutine particles_allocate use m_io implicit none write(out,*) 'Particle allocation' call flush(out) ! allocate the arrays if (.not.allocated(xyzp)) then allocate(xyzp(3,nptot),uvwp(3,nptot),ipart(nptot),itmp_part(nptot),& myid_part(nptot),& xyzp1(3,nptot),uvwp1(3,nptot),ipart1(nptot),stat=ierr) if (ierr.ne.0) stop 'particle allocation' else write(out,*) 'arrays xyzp etc are already allocated.' call flush(out) end if ! allocating temp arrays for CINT interpolation if (particles_tracking_scheme.eq.2) then write(out,*) 'Allocating temp arrays for CINT interpolation' call flush(out) if (.not.allocated(wrk1p)) then allocate(wrk1p(nx,ny), wrk2p(nx,ny), wrk3p(nx,ny),stat=ierr) wrk1p = zip; wrk2p = zip; wrk3p = zip; if (ierr.ne.0) stop 'allocation of particle temp arrays' else write(out,*) 'arrays for CINT are already allocated.' call flush(out) end if end if return end subroutine particles_allocate !================================================================================ !================================================================================ !================================================================================ subroutine particles_deallocate implicit none write(out,*) 'Particle deallocation.' call flush(out) ! allocate the arrays if (allocated(xyzp)) deallocate(xyzp,uvwp,ipart,itmp_part,myid_part) if (allocated(wrk1p)) deallocate(wrk1p,wrk2p,wrk3p,xyzp1,uvwp1,ipart1) return end subroutine particles_deallocate !================================================================================ !================================================================================ !================================================================================ ! PARTICLES_GET_MYID - fills the array myid_part, which contains the IDs of ! the processes that house the particles. !================================================================================ subroutine particles_get_myid implicit none integer :: n ! filling up the array which contains id's of processes that are ! responsible for each particle itmp_part = 0 myid_part = 0 if (np.gt.0) then do n = 1,np itmp_part( ipart(n) ) = myid end do end if call MPI_REDUCE(itmp_part,myid_part,nptot,MPI_INTEGER4,MPI_SUM,master,MPI_COMM_TASK,ierr) call MPI_BCAST(myid_part,nptot,MPI_INTEGER4,master,MPI_COMM_TASK,ierr) return end subroutine particles_get_myid !================================================================================ !================================================================================ !================================================================================ ! PARTICLES_INIT - initializes the particles !================================================================================ subroutine particles_init !-------------------------------------------------------------------------------- ! particle initialization: ! - aray allocation ! - reading coordinates of particles from the file .pt ! - shifing the particles so hey are in [0:2*pi) ! - distributing the particles among CPUs !-------------------------------------------------------------------------------- use m_filter_xfftw implicit none logical :: init_start = .false., init_internally = .false. integer :: i,j,k,npthird,n real*8 :: sctmp !-------------------------------------------------------------------------------- ! first, some foolproofing int_particles = .false. ! if the total # of particles is zero, return without initialization if (nptot.eq.0) then write(out,*) ' nptot = 0, no particles initialized.' call flush(out) return end if !-------------------------------------------------------------------------------- ! choose which way to initialize particles: ! 1. read from file .pt (manually given coordinates) ! 2. generate uniformly spaced particles in the computational domain ! (in this case the number of particles is taken to be nptot=2^n, where ! n is such that nptot is close to what is specified in the input file) !-------------------------------------------------------------------------------- ! if the file pt. is there, then read particles from it and return. ! whatever is specified in the pt., supercedes anything else. inquire(file='pt.'//file_ext,exist=there) if (there) then write(out,*) 'Detected file pt.'//file_ext write(out,*) 'Particles will be read from the file.' call particles_restart_read_binary return end if ! if the restart file is not there, we assume that it's a ! brand new particle simulation so we read the particle ! coordinates from the "run_name".pt file ! --------- OR ---------- ! define particles as uiformly distributed over the domain ! check if the .pt is in the directory inquire(file=run_name//'.pt',exist=there) if (there) then write(out,*) 'Reading particles from the file '//run_name//'.pt' call flush(out) init_start = .true. else ! if not reading from the file, redefine nptot init_internally = .true. npthird = int(dble(nptot)**(1.d0/3.d0)) if (nptot-npthird**3.ge.(npthird+1)**3-nptot) npthird = npthird + 1 nptot = npthird**3 nptot = min(nptot,2**30) write(out,*) 'Redefined nptot to be a perfect cube:',nptot if (nptot.eq.2**30) write(out,*) "REACHED MAXIMUM: 2**30" call flush(out) end if write(out,*) 'Initializing',nptot,' particles' call flush(out) call particles_allocate if (myid.eq.0) then ! read the whole array of particles if (init_start) then open(99,file=run_name//'.pt') read(99,*,ERR=9000,END=9000) ((xyzp(i,j),i=1,3),j=1,nptot) close(99) end if ! ... or define it manually sctmp = two*PI / real(npthird,8) if (init_internally) then n = 0 do k = 1,npthird do j = 1,npthird do i = 1,npthird n = n + 1 xyzp(1,n) = (dble(i-1)+half) * sctmp xyzp(2,n) = (dble(j-1)+half) * sctmp xyzp(3,n) = (dble(k-1)+half) * sctmp end do end do end do end if end if ! Broadcast the array of particles call MPI_BCAST(xyzp,nptot*3,MPI_REAL8,0,MPI_COMM_TASK,mpi_err) ! definition of particles' identifiers do i = 1,nptot ipart(i) = i end do !!$ write(out,*) 'Particles read in:' !!$ write(out,'(3f20.10)') ((xyzp(i,j),i=1,3),j=1,nptot) !!$ call flush(out) ! REARRANGING THE PARTICLES ! makign sure the xyz are in [0,2*pi) do j=1,nptot do i = 1,3 if (xyzp(i,j) .ge. 2.0d0*PI .or. xyzp(i,j) .lt. 0.0d0) then xyzp(i,j) = xyzp(i,j) - 2.0d0*PI * dble(floor(xyzp(i,j)/(2.0d0*PI))) end if end do end do !!$ write(out,*) 'Particles normalized:',nptot !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,nptot) ! removing particles that do not correspond to this slab np = nptot j = 1 zmin = dble(myid*nz) * 2.0d0*PI / dble(nz_all) zmax = dble((myid+1)*nz) * 2.0d0*PI / dble(nz_all) do while (j.le.np) if (xyzp(3,j).ge.zmin .and. xyzp(3,j).lt.zmax) then j = j + 1 else xyzp(:,j) = xyzp(:,np) ipart(j) = ipart(np) np = np-1 end if end do !!$ write(out,*) 'Particles left in:',np !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,np) ! Recalculating the particles' coordinates in grid cells do j = 1,np xyzp(1,j) = xyzp(1,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 xyzp(2,j) = xyzp(2,j) / (2.0d0 * PI / dble(ny)) + 1.0d0 xyzp(3,j) = xyzp(3,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 - dble(myid*nz) end do write(out,*) 'Particles left in, rescaled:',np !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,np) ! initializing uvwp with -999.00 uvwp = -999.0d0 !!$ ! making a subdirectory "particles" !!$ if (iammaster) call system('mkdir -p particles') ! in case if the particles are advected by locally averaged field ! (particles_filter_size > 0), define the filter if (particles_filter_size > 0.d0) call filter_xfftw_init return 9000 write(out,*) '*** PARTICLES_INIT: ERROR in reading file .pt !!!' stop end subroutine particles_init !=========================================================================== !=========================================================================== ! Subroutine that writes out the particles coordinates in a file ! during the restart write !=========================================================================== subroutine particles_restart_write implicit none integer(kind=MPI_INTEGER_KIND) :: id, count, my_status(MPI_STATUS_SIZE) integer*4 :: np_rec integer :: np_cur, ntmp if (.not.int_particles) return write(out,*) 'writing the particle restart file' call flush(out) ! ------------------------------------------------------------------------------ ! the slave nodes send their particles' IDs and coordinates to the master node ! the master node assembles the particles IDs and coordinates in one array and ! writes them out ! ------------------------------------------------------------------------------ if (myid.ne.0) then ! First send the number of particles that the process holds to the master process np_rec = np call MPI_SEND(np_rec,1,MPI_INTEGER4,0,3*myid,MPI_COMM_TASK,mpi_err) ! write(out,*) 'Sent np=',np_rec,' to master',mpi_err ! call flush(out) ! if the number of particles is positive, send their IDs and coordinates if (np.gt.0) then count = np call MPI_SEND(ipart(1),count,MPI_INTEGER4,0,3*myid+1,MPI_COMM_TASK,mpi_err) ! write(out,*) 'Sent ipart to master',mpi_err ! call flush(out) count = 3*np call MPI_SEND(xyzp(1,1),count,MPI_REAL8,0,3*myid+2,MPI_COMM_TASK,mpi_err) ! write(out,*) 'Sent xyzp to master',mpi_err ! call flush(out) end if else np_cur = np+1 do id = 1,numprocs-1 ! master node receives the number of particles from slave node call MPI_RECV(np_rec,1,MPI_INTEGER4,id,3*id,MPI_COMM_TASK,my_status,mpi_err) ! write(out,*) 'Received ',np_rec,' from ',id,':',mpi_err ! call flush(out) ! if there are any particles to receive, receive the IDs and coordinates if (np_rec.gt.0) then count = np_rec call MPI_RECV(ipart(np_cur),count,MPI_INTEGER4,id,3*id+1,MPI_COMM_TASK,my_status,mpi_err) ! write(out,*) 'Received ipart from ',id,':',mpi_err ! call flush(out) count = 3 * np_rec call MPI_RECV(xyzp(1,np_cur),count,MPI_REAL8,id,3*id+2,MPI_COMM_TASK,my_status,mpi_err) ! write(out,*) 'Received xyzp from ',id,':',mpi_err ! call flush(out) ! since the xyzp contains particle coordinates in terms of the cells, the third ! coordinate should be augmented by id*nz. this is done so later we can just ! multiply the third coordinate by dz to get the real particle coordinate do ntmp = np_cur,np_cur+np_rec xyzp(3,ntmp) = xyzp(3,ntmp) + dble(nz*id) end do end if np_cur = np_cur + np_rec end do end if ! ------------------------------------------------------------------------------ ! Now the master node opens the particle restart file .pt. ! and writes the particles with their IDs in it. ! ------------------------------------------------------------------------------ if (myid.eq.0) then open(89,file=run_name//'.pt.'//file_ext) do ntmp = 1,nptot write(89,'(i6,x,10d16.8)') ipart(ntmp),& (xyzp(1,ntmp)-1.d0)*dx,(xyzp(2,ntmp)-1.d0)*dy,(xyzp(3,ntmp)-1.d0)*dz end do close(89,status='keep') write(out,*) 'particle restart file written.' call flush(out) end if return end subroutine particles_restart_write !=========================================================================== !=========================================================================== ! Subroutine that reads in the particles from the particle restart file !=========================================================================== subroutine particles_restart_read implicit none logical :: there character*80 :: fname integer :: n,i,j if (nptot.eq.0) return ! allocate the arrays if (.not.allocated(xyzp)) then allocate(xyzp(3,nptot),uvwp(3,nptot),ipart(nptot),itmp_part(nptot),& myid_part(nptot),stat=ierr) if (ierr.ne.0) stop 'particle allocation' else write(out,*) 'PARTICLES_RESTART_READ: xyzp allocated already!' call flush(out) end if fname = run_name//'.pt.'//file_ext write(out,*) 'Reading the particles from file '//fname call flush(out) !------------------------------------------------------------------------- ! Reading and broadcasting the particles' IDs and coordinates !------------------------------------------------------------------------- inquire(file=fname,exist=there) if(.not.there) then write(out,*) 'could not find the particle restart file ',fname write(out,*) 'exiting.' call flush(out) call my_exit(-1) end if if (myid.eq.0) then open(81,file=fname) do n = 1,nptot read(81,*,ERR=9000,END=9000) ipart(n),xyzp(1:3,n) end do close(81) end if call MPI_BCAST(ipart,nptot,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err) call MPI_BCAST(xyzp,nptot*3,MPI_REAL8,0,MPI_COMM_TASK,mpi_err) write(out,*) 'Particles normalized:',nptot do j = 1,nptot write(out,'(i5,3f20.10)') ipart(j),(xyzp(i,j),i=1,3) end do !------------------------------------------------------------------------- ! removing particles that do not correspond to this slab !------------------------------------------------------------------------- write(out,*) 'removing particles that do not correspond to this slab' call flush(out) np = nptot j = 1 zmin = dble(myid*nz) * 2.0d0*PI / dble(nz_all) zmax = dble((myid+1)*nz) * 2.0d0*PI / dble(nz_all) do while (j.le.np) if (xyzp(3,j).ge.zmin .and. xyzp(3,j).lt.zmax) then j = j + 1 else xyzp(:,j) = xyzp(:,np) ipart(j) = ipart(np) np = np-1 end if end do write(out,*) 'Particles left in:',np write(out,*) 'Particles normalized:',nptot do j = 1,nptot write(out,'(i5,3f20.10)') ipart(j),(xyzp(i,j),i=1,3) end do ! Recalculating the particles' coordinates in grid cells do j = 1,np xyzp(1,j) = xyzp(1,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 xyzp(2,j) = xyzp(2,j) / (2.0d0 * PI / dble(ny)) + 1.0d0 xyzp(3,j) = xyzp(3,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 - dble(myid*nz) end do write(out,*) 'Particles left in, rescaled:',np write(out,*) 'Particles normalized:',nptot do j = 1,nptot write(out,'(i5,3f20.10)') ipart(j),(xyzp(i,j),i=1,3) end do ! initializing uvwp with -999.00 uvwp = -999.0d0 return 9000 write(out,*) 'Error reading particle restart file ',fname stop end subroutine particles_restart_read !================================================================================ !================================================================================ !================================================================================ ! Subroutine that writes out the particles coordinates in a BINARY restart file !================================================================================ subroutine particles_restart_write_binary implicit none integer(kind=MPI_INTEGER_KIND) :: id, count, my_status(MPI_STATUS_SIZE) integer*4 :: np_rec integer :: np_cur, ntmp ! if the last iteration at which the particles were written out is the ! current iteration, then skip the writing if (particles_last_dump .eq. itime) return ! otherwise proceed and redefine the particles_last_dump particles_last_dump = itime !!$ write(out,*) 'writing the BINARY particle restart file' !!$ call flush(out) ! ------------------------------------------------------------------------------ ! the slave nodes send their particles' IDs and coordinates to the master node ! the master node assembles the particles IDs and coordinates in one array and ! writes them out ! ------------------------------------------------------------------------------ if (myid.ne.0) then ! First send the number of particles that the process holds to the master process np_rec = np call MPI_SEND(np_rec,1,MPI_INTEGER4,0,3*myid,MPI_COMM_TASK,mpi_err) ! if the number of particles is positive, send their IDs and coordinates if (np.gt.0) then count = np call MPI_SEND(ipart(1),count,MPI_INTEGER4,0,3*myid+1,MPI_COMM_TASK,mpi_err) count = 3*np call MPI_SEND(xyzp(1,1),count,MPI_REAL8,0,3*myid+2,MPI_COMM_TASK,mpi_err) end if else np_cur = np+1 do id = 1,numprocs-1 ! master node receives the number of particles from slave node call MPI_RECV(np_rec,1,MPI_INTEGER4,id,3*id,MPI_COMM_TASK,my_status,mpi_err) ! if there are any particles to receive, receive the IDs and coordinates if (np_rec.gt.0) then count = np_rec call MPI_RECV(ipart(np_cur),count,MPI_INTEGER4,id,3*id+1,MPI_COMM_TASK,my_status,mpi_err) count = 3 * np_rec call MPI_RECV(xyzp(1,np_cur),count,MPI_REAL8,id,3*id+2,MPI_COMM_TASK,my_status,mpi_err) ! since the xyzp contains particle coordinates in terms of the cells, the third ! coordinate should be augmented by id*nz. this is done so later we can just ! multiply the third coordinate by dz to get the real particle coordinate do ntmp = np_cur,np_cur+np_rec xyzp(3,ntmp) = xyzp(3,ntmp) + dble(nz*id) end do end if np_cur = np_cur + np_rec end do end if ! ------------------------------------------------------------------------------ ! Now the master node opens the particle BINARY file pt. ! and writes the particles with their IDs in it. ! ------------------------------------------------------------------------------ if (myid.eq.0) then !! open(89,file='pt.'//file_ext,form='binary') open(89,file='pt.'//file_ext,form='unformatted',access='stream') write(89) nptot do ntmp = 1,nptot write(89) ipart(ntmp),& (xyzp(1,ntmp)-1.d0)*dx,(xyzp(2,ntmp)-1.d0)*dy,(xyzp(3,ntmp)-1.d0)*dz end do ! writing time stamp in the particle output file write(89) TIME close(89,status='keep') write(out,'("particle BINARY file written:",i7,e15.6)') itime, time call flush(out) end if return end subroutine particles_restart_write_binary !================================================================================ !=========================================================================== ! Subroutine that reads in the particles from the particle restart file !=========================================================================== subroutine particles_restart_read_binary implicit none integer(kind=MPI_INTEGER_KIND) :: nptot1 logical :: there character*80 :: fname integer :: n,i,j fname = 'pt.'//file_ext write(out,*) 'Reading the particles from binary file '//fname call flush(out) !------------------------------------------------------------------------- ! Reading and broadcasting the particles' IDs and coordinates !------------------------------------------------------------------------- ! Checking that the number of particles in the binary file is correct if (myid.eq.0) then open(81,file=fname,form='unformatted',access='stream') read(81) nptot end if call MPI_BCAST(nptot,1,MPI_INTEGER,0,MPI_COMM_TASK,mpi_err) write(out,*) 'nptot = ',nptot call flush(out) ! allocating arrays call particles_allocate ! Readign and broadcasting particles if (myid.eq.0) then do n = 1,nptot read(81,ERR=9000,END=9000) ipart(n),xyzp(1:3,n) end do close(81) end if call MPI_BCAST(ipart,nptot,MPI_INTEGER4,0,MPI_COMM_TASK,mpi_err) call MPI_BCAST(xyzp,nptot*3,MPI_REAL8,0,MPI_COMM_TASK,mpi_err) write(out,*) 'Particles normalized:',nptot !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,nptot) !------------------------------------------------------------------------- ! removing particles that do not correspond to this slab !------------------------------------------------------------------------- write(out,*) 'removing particles that do not correspond to this slab' call flush(out) np = nptot j = 1 zmin = dble(myid*nz) * 2.0d0*PI / dble(nz_all) zmax = dble((myid+1)*nz) * 2.0d0*PI / dble(nz_all) do while (j.le.np) if (xyzp(3,j).ge.zmin .and. xyzp(3,j).lt.zmax) then j = j + 1 else xyzp(:,j) = xyzp(:,np) ipart(j) = ipart(np) np = np-1 end if end do write(out,*) 'Particles left in:',np !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,np) ! Recalculating the particles' coordinates in grid cells do j = 1,np xyzp(1,j) = xyzp(1,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 xyzp(2,j) = xyzp(2,j) / (2.0d0 * PI / dble(ny)) + 1.0d0 xyzp(3,j) = xyzp(3,j) / (2.0d0 * PI / dble(nx)) + 1.0d0 - dble(myid*nz) end do !!$ write(out,*) 'Particles left in, rescaled:',np !!$ write(out,'(i5,3f20.10)') ((ipart(j),(xyzp(i,j),i=1,3)),j=1,np) ! initializing uvwp with -999.00 uvwp = -999.0d0 return 9000 write(out,*) 'Error reading particle restart file ',fname stop end subroutine particles_restart_read_binary !================================================================================ !================================================================================ !================================================================================ subroutine particles_move use m_fields use m_work use x_fftw use m_filter_xfftw implicit none integer :: m, n if (task.ne.'parts') return ! if the particles are advected by locally averaged velocities ! instead of fully resolved velocities, then we need to filter ! (locally average) these velocities. The indicator of this case is ! particles_filter_size > 0.d0 ! in this case the velocities that are sent to the "parts" part of the code ! are given in the Fourier space. So we need to ! 1) filter them ! 2) transform them back to real space locally_averaged_velocity: if (particles_filter_size .gt. 0.d0) then ! filter each field and transform to real space do n = 1,3 call filter_xfftw_fields(n) call xFFT3d_fields(-1,n) end do end if locally_averaged_velocity ! now shifting particles according to the interpolation scheme and advection scheme select case (particles_tracking_scheme) case (0) call particles_interpolate_trilinear case (1) call particles_interpolate_cint case default stop 'wrong particles_tracking_scheme' end select call particles_update_slabs return end subroutine particles_move !================================================================================ ! PARTICLE VELOCITY INTERPOLATION AND MOVING ! ! This is set up in the following way: ! ! The interpolation subroutines do the following: ! 1) find particles' velocities using some interpolation methof ! 2) move the particles according to the numerical scheme ! (Euler or Adams-Bashforth) ! 3) Calculate the number of particles that went out of the slab: ! np_send_d - # of particles that needs to be sent "down", to slab # myid-1 ! np_send_u - # of particles that needs to be sent "up", to slab # myid+1 ! ! Then the subroutine particles_update_slabs should be called, which ! 1) rearranges particles according to their destination (here, down, up) ! 2) inform the neighbor slabs about the number of particles to send ! 3) transfer particles between the slabs ! ! May the Force be with us. !================================================================================ !================================================================================ subroutine particles_interpolate_trilinear use m_fields use m_work implicit none integer :: i,j,k,i1,j1,k1,n real*8 :: xp,yp,zp,up,vp,wp real*8 :: c1,c2,c3,c4,c5,c6,c7,c8 ! note that we assume that fields1...3 contain REAL velocities (in X-space) ! to interpolate velocities of the particles that are on the upper z-edge of ! the slice, we need the first layer of velocities from the neibour process. ! This is done by sending them to the process with number (myid-1) and storing ! in the first three slices (k=1,3) of wrk0 id_to = mod(myid-1+numprocs,numprocs) id_from = mod(myid+1,numprocs) ! sending u sendtag = 3 * myid recvtag = 3 * id_from count = (nx+2)*ny call MPI_IRECV( wrk(1,1,1,0), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,1), count, MPI_REAL8, id_to, sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) ! sending v sendtag = 3 * myid + 1 recvtag = 3 * id_from + 1 count = (nx+2)*ny call MPI_IRECV( wrk(1,1,2,0), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,2), count, MPI_REAL8, id_to, sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) ! sending w sendtag = 3 * myid + 2 recvtag = 3 * id_from + 2 count = (nx+2)*ny call MPI_IRECV( wrk(1,1,3,0), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,3), count, MPI_REAL8, id_to, sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) ! # of particles to be sent/received to/from left and right neibors np_send_d = 0 np_send_u = 0 np_get_d = 0 np_get_u = 0 ! interpolating the velocities and moving do n = 1,np ! getting the cell where the particle is at i = floor(xyzp(1,n)) j = floor(xyzp(2,n)) k = floor(xyzp(3,n)) ! getting fractional coordinates of the particle xp = xyzp(1,n) - dble(i) yp = xyzp(2,n) - dble(j) zp = xyzp(3,n) - dble(k) ! cpefficients for trilinear interpolation c1 = (one-xp) * (one-yp) * (one-zp) c2 = (xp) * (one-yp) * (one-zp) c3 = (one-xp) * (yp) * (one-zp) c4 = (xp) * (yp) * (one-zp) c5 = (one-xp) * (one-yp) * (zp) c6 = (xp) * (one-yp) * (zp) c7 = (one-xp) * (yp) * (zp) c8 = (xp) * (yp) * (zp) ! accounting for periodicity in x and y directions if (i.eq.nx) then i1 = 1 else i1 = i + 1 end if if (j.eq.ny) then j1 = 1 else j1 = j + 1 end if ! interpolating velocities to get the velocity of the particle up = c1 * fields(i,j ,k,1) + c2 * fields(i1,j ,k,1) + & c3 * fields(i,j1,k,1) + c4 * fields(i1,j1,k,1) vp = c1 * fields(i,j ,k,2) + c2 * fields(i1,j ,k,2) + & c3 * fields(i,j1,k,2) + c4 * fields(i1,j1,k,2) wp = c1 * fields(i,j ,k,3) + c2 * fields(i1,j ,k,3) + & c3 * fields(i,j1,k,3) + c4 * fields(i1,j1,k,3) if (k.eq.nz) then ! if the particle is in the last layer, add velocities from wrk0 up = up + & c5 * wrk(i,j ,1,0) + c6 * wrk(i1,j ,1,0) + & c7 * wrk(i,j1,1,0) + c8 * wrk(i1,j1,1,0) vp = vp + & c5 * wrk(i,j ,2,0) + c6 * wrk(i1,j ,2,0) + & c7 * wrk(i,j1,2,0) + c8 * wrk(i1,j1,2,0) wp = wp + & c5 * wrk(i,j ,3,0) + c6 * wrk(i1,j ,3,0) + & c7 * wrk(i,j1,3,0) + c8 * wrk(i1,j1,3,0) else ! if the particle is inside the slice, straightforward interpolation up = up + & c5 * fields(i,j ,k+1,1) + c6 * fields(i1,j ,k+1,1) + & c7 * fields(i,j1,k+1,1) + c8 * fields(i1,j1,k+1,1) vp = vp + & c5 * fields(i,j ,k+1,2) + c6 * fields(i1,j ,k+1,2) + & c7 * fields(i,j1,k+1,2) + c8 * fields(i1,j1,k+1,2) wp = wp + & c5 * fields(i,j ,k+1,3) + c6 * fields(i1,j ,k+1,3) + & c7 * fields(i,j1,k+1,3) + c8 * fields(i1,j1,k+1,3) end if ! now move the particle if (uvwp(1,n).eq.-999.0d0) then ! if the second slice is unavailable, use Euler xp = xp + dt * up / dx yp = yp + dt * vp / dy zp = zp + dt * wp / dz else ! if there is the second slice, use Adams-Bashforth xp = xp + (1.50d0*up - 0.50d0*uvwp(1,n)) * dt/dx yp = yp + (1.50d0*vp - 0.50d0*uvwp(2,n)) * dt/dy zp = zp + (1.50d0*wp - 0.50d0*uvwp(3,n)) * dt/dz end if ! update the velocities uvwp(1,n) = up uvwp(2,n) = vp uvwp(3,n) = wp ! apply periodicity in x- and y-directions ! also count the # of particles to be sent to neibour slices ! note that we assume that the particle cannot move further than one cell ! from its curent position ! x direction if (xp.gt.1.0d0) then xp = xp - 1.0d0 i = i+1 if (i.gt.nx) i = i - nx end if if (xp.lt.0.0d0) then xp = xp + 1.0d0 i = i-1 if (i.lt.1) i = i + nx end if ! y direction if (yp.gt.1.0d0) then yp = yp - 1.0d0 j = j+1 if (j.gt.ny) j = j - ny end if if (yp.lt.0.0d0) then yp = yp + 1.0d0 j = j-1 if (j.lt.1) j = j + ny end if ! z direction if (zp.gt.1.0d0) then zp = zp - 1.0d0 k = k+1 if (k.gt.nz) np_send_u = np_send_u + 1 end if if (zp.lt.0.0d0) then zp = zp + 1.0d0 k = k-1 if (k.lt.1) np_send_d = np_send_d + 1 end if ! updating the particle coordinates. ! note that z-coordinate can be outside of the slab xyzp(1,n) = xp + dble(i) xyzp(2,n) = yp + dble(j) xyzp(3,n) = zp + dble(k) end do return end subroutine particles_interpolate_trilinear !================================================================================ !================================================================================ ! Moving Lagrangian particles. The velocity is interpolated using tricubic ! interpolation method, found in Wikipedia :) !================================================================================ subroutine particles_interpolate_cint use m_fields use m_work implicit none ! local variables integer :: ip,jp,kp,i,j,k,i1,j1,k1,n real*8 :: xp,yp,zp,up,vp,wp real*8 :: sctmp real*8 :: uloc(4,4,4), vloc(4,4,4), wloc(4,4,4) ! real*8 :: cint_3d real*8 :: r ! making sure that the thickness of the slabs is mroe than 1 slice if (nz.lt.2) then write(out,*) "PARTICLES_MOVE_CINT: nz is less than 2, stopping" call flush(out) stop 'particles_move_cint *** nz less than 2' end if ! assuming that array fields1...3 has the velocities in X-space ! because this method involves the adjacent 26 cells (a 4x4 cube of points ! in which the particle occupies the central cell), we need to pass the ! first two slices "down" and the last single slice (nz) "up" ! sending/receiving the single slice ! sending "up" ! receiving from "down" count = (nx+2) * ny id_from = mod(myid - 1 + numprocs, numprocs); id_to = mod(myid + 1 , numprocs); sendtag = 6 * myid + 0 recvtag = 6 * id_from + 0 call MPI_IRECV( wrk(1,1,1, 1), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,nz,1), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) sendtag = 6 * myid + 1 recvtag = 6 * id_from + 1 call MPI_IRECV( wrk(1,1,1, 2), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,nz,2), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) sendtag = 6 * myid + 2 recvtag = 6 * id_from + 2 call MPI_IRECV( wrk(1,1,1, 3), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,nz,3), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) ! sending/receiving two slicee ! sending "down" ! receiving from "up" count = 2 * (nx+2) * ny id_from = mod(myid + 1 , numprocs); id_to = mod(myid - 1 + numprocs, numprocs); sendtag = 6 * myid + 3 recvtag = 6 * id_from + 3 call MPI_IRECV( wrk(1,1,2,1), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,1), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) sendtag = 6 * myid + 4 recvtag = 6 * id_from + 4 call MPI_IRECV( wrk(1,1,2,2), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,2), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) sendtag = 6 * myid + 5 recvtag = 6 * id_from + 5 call MPI_IRECV( wrk(1,1,2,3), count, MPI_REAL8, id_from, recvtag, MPI_COMM_TASK, request, mpi_err) call MPI_SEND(fields(1,1,1,3), count, MPI_REAL8, id_to , sendtag, MPI_COMM_TASK, mpi_err) call MPI_WAIT(request, mpi_status, mpi_err) ! # of particles to be sent/received to/from "down" and "up" np_send_u = 0 np_send_d = 0 np_get_u = 0 np_get_d = 0 do n = 1,np ! getting the cell where the particle is at ip = floor(xyzp(1,n)) jp = floor(xyzp(2,n)) kp = floor(xyzp(3,n)) ! getting fractional coordinates of the particle xp = xyzp(1,n) - dble(ip) yp = xyzp(2,n) - dble(jp) zp = xyzp(3,n) - dble(kp) ! ------------------------------------------------------------------- ! Getting the velocities of the particles by tricubic interpolation ! (for particularities, see notes from 12/13/07) ! fillin out arrays uloc, vloc, and wloc(4,4,4) for tricubic interpolation ! taking into account periodicity in x and y. do k = 1,4 k1 = kp - 2 + k do j = 1,4 j1 = jp - 2 + j if (j1.lt.1 ) j1 = j1 + ny if (j1.gt.ny) j1 = j1 - ny do i = 1,4 i1 = ip - 2 + i if (i1.lt.1 ) i1 = i1 + nx if (i1.gt.nx) i1 = i1 - nx if (k1.lt.1) then uloc(i,j,k) = wrk(i1,j1,1,1) vloc(i,j,k) = wrk(i1,j1,1,2) wloc(i,j,k) = wrk(i1,j1,1,3) else if (k1.gt.nz) then uloc(i,j,k) = wrk(i1,j1,k1-nz+1,1) vloc(i,j,k) = wrk(i1,j1,k1-nz+1,2) wloc(i,j,k) = wrk(i1,j1,k1-nz+1,3) else uloc(i,j,k) = fields(i1,j1,k1,1) vloc(i,j,k) = fields(i1,j1,k1,2) wloc(i,j,k) = fields(i1,j1,k1,3) end if end do end do end do ! now we have arrays uloc, vloc, wloc and can interpolate up = cint_3d(uloc,xp,yp,zp) vp = cint_3d(vloc,xp,yp,zp) wp = cint_3d(wloc,xp,yp,zp) ! now move the particle if (uvwp(1,n).eq.-999.0d0) then ! if the second slice is unavailable, use Euler xp = xp + dt * up / dx yp = yp + dt * vp / dy zp = zp + dt * wp / dz else ! if there is the second slice, use Adams-Bashforth xp = xp + (1.50d0*up - 0.50d0*uvwp(1,n)) * dt/dx yp = yp + (1.50d0*vp - 0.50d0*uvwp(2,n)) * dt/dy zp = zp + (1.50d0*wp - 0.50d0*uvwp(3,n)) * dt/dz end if uvwp(1,n) = up uvwp(2,n) = vp uvwp(3,n) = wp ! apply periodicity in x- and y-directions ! also count the # of particles to be sent to neibour slices ! note that we assume that the particle cannot move further than one cell ! from its curent position ! x direction if (xp.gt.1.0d0) then xp = xp - 1.0d0 ip = ip+1 if (ip.gt.nx) ip = ip - nx end if if (xp.lt.0.0d0) then xp = xp + 1.0d0 ip = ip-1 if (ip.lt.1) ip = ip + nx end if ! y direction if (yp.gt.1.0d0) then yp = yp - 1.0d0 jp = jp+1 if (jp.gt.ny) jp = jp - ny end if if (yp.lt.0.0d0) then yp = yp + 1.0d0 jp = jp-1 if (jp.lt.1) jp = jp + ny end if ! z direction if (zp.gt.1.0d0) then zp = zp - 1.0d0 kp = kp+1 if (kp.gt.nz) np_send_u = np_send_u + 1 end if if (zp.lt.0.0d0) then zp = zp + 1.0d0 kp = kp-1 if (kp.lt.1) np_send_d = np_send_d + 1 end if ! updating the particle coordinates. ! note that z-coordinate can be outside of the slice xyzp(1,n) = xp + dble(ip) xyzp(2,n) = yp + dble(jp) xyzp(3,n) = zp + dble(kp) end do return end subroutine particles_interpolate_cint !================================================================================ !================================================================================ subroutine particles_update_slabs use m_openmpi use m_work implicit none integer (kind=MPI_INTEGER_KIND) :: id_down, id_up integer :: i, j, n, iii id_down = mod(myid-1+numprocs,numprocs) id_up = mod(myid+1,numprocs) !!$ write(out,*) "np_sends:", np_send_d, np_send_u !!$ call flush(out) ! let the neighbors know how many particles went to them call MPI_SEND(np_send_d,1,MPI_INTEGER4,id_down,2*myid ,MPI_COMM_TASK,mpi_err) call MPI_SEND(np_send_u,1,MPI_INTEGER4,id_up ,2*myid+1,MPI_COMM_TASK,mpi_err) ! learn how many particles will be transferred to this slice from neighbors call MPI_RECV(np_get_d,1,MPI_INTEGER4,id_down,2*id_down+1,MPI_COMM_TASK,mpi_status,mpi_err) call MPI_RECV(np_get_u,1,MPI_INTEGER4,id_up ,2*id_up ,MPI_COMM_TASK,mpi_status,mpi_err) !!$ write(out,*) "np_gets:", np_get_d, np_get_u !!$ call flush(out) ! rearranging particles in the following order: ! 1. particles that stay ! 2. particles that leave down ! 3. particles that leave up if (np_send_u.gt.0) then i = 1 n = np do while (i.le.n) if (xyzp(3,i).ge.dble(nz+1)) then wrk(1:3,1,1,0) = xyzp(1:3,i) wrk(4:6,1,1,0) = uvwp(1:3,i) j = ipart(i) xyzp(:,i) = xyzp(:,n) uvwp(:,i) = uvwp(:,n) ipart(i) = ipart(n) xyzp(1:3,n) = wrk(1:3,1,1,0) uvwp(1:3,n) = wrk(4:6,1,1,0) xyzp(3,n) = xyzp(3,n) - dble(nz) ipart(n) = j n = n - 1 else i = i + 1 end if end do end if if (np_send_d.gt.0) then i = 1 n = np - np_send_u do while (i.le.n) if (xyzp(3,i).lt.1.0d0) then wrk(1:3,1,1,0) = xyzp(1:3,i) wrk(4:6,1,1,0) = uvwp(1:3,i) j = ipart(i) xyzp(:,i) = xyzp(:,n) uvwp(:,i) = uvwp(:,n) ipart(i) = ipart(n) xyzp(1:3,n) = wrk(1:3,1,1,0) uvwp(1:3,n) = wrk(4:6,1,1,0) xyzp(3,n) = xyzp(3,n) + dble(nz) ipart(n) = j n = n - 1 else i = i + 1 end if end do end if ! sending and receiving particiles from adjacent slabs ! extra arrays (a dumb way to avoid deadlocking in send/receive when ! number of particles is large and the send/receive is not completed in the ! single handshake between the two processes) xyzp1 = xyzp uvwp1 = uvwp ipart1 = ipart np1 = np np = np - np_send_d - np_send_u ! posting receive from down if (np_get_d.gt.0) then call MPI_IRECV(xyzp(1,np+1),3*np_get_d,MPI_REAL8 ,id_down,6*id_down+0,MPI_COMM_TASK,request1,mpi_err) call MPI_IRECV(uvwp(1,np+1),3*np_get_d,MPI_REAL8 ,id_down,6*id_down+1,MPI_COMM_TASK,request2,mpi_err) call MPI_IRECV( ipart(np+1),np_get_d,MPI_INTEGER8,id_down,6*id_down+2,MPI_COMM_TASK,request3,mpi_err) !!$ write(out,*) 'posted receive from ',id_down !!$ call flush(out) end if ! sending particles to up if (np_send_u.gt.0) then !!$ write(out,*) 'xyzp1' !!$ write(out,"(3e15.6)") (xyzp1(:,iii),iii=1,np) !!$ write(out,*) 'uvwp1' !!$ write(out,"(3e15.6)") (uvwp1(:,iii),iii=1,np) !!$ write(out,*) 'ipart1' !!$ write(out,"(i5)") (ipart(iii),iii=1,np) !!$ call flush(out) call MPI_SEND(xyzp1(1,np1-np_send_u+1),3*np_send_u,MPI_REAL8 ,id_up,6*myid+0,MPI_COMM_TASK,mpi_err) call MPI_SEND(uvwp1(1,np1-np_send_u+1),3*np_send_u,MPI_REAL8 ,id_up,6*myid+1,MPI_COMM_TASK,mpi_err) call MPI_SEND( ipart1(np1-np_send_u+1),np_send_u,MPI_INTEGER4,id_up,6*myid+2,MPI_COMM_TASK,mpi_err) !!$ write(out,*) 'sent',np_send_u,' to the right:', ipart1(np1-np_send_u+1:np1) !!$ call flush(out) np1 = np1 - np_send_u end if ! completing the receive from the down if (np_get_d.gt.0) then call MPI_WAIT(request1,mpi_status,mpi_err) call MPI_WAIT(request2,mpi_status,mpi_err) call MPI_WAIT(request3,mpi_status,mpi_err) !!$ write(out,*) 'received',np_get_d,' from the left:', ipart(np+1:np+np_get_d) !!$ call flush(out) np = np + np_get_d end if ! posting receive from the up if (np_get_u.gt.0) then call MPI_IRECV(xyzp(1,np+1),3*np_get_u,MPI_REAL8 ,id_up,6*id_up+3,MPI_COMM_TASK,request1,mpi_err) call MPI_IRECV(uvwp(1,np+1),3*np_get_u,MPI_REAL8 ,id_up,6*id_up+4,MPI_COMM_TASK,request2,mpi_err) call MPI_IRECV( ipart(np+1),np_get_u ,MPI_INTEGER8,id_up,6*id_up+5,MPI_COMM_TASK,request3,mpi_err) !!$ write(out,*) 'posted receive from ',id_up !!$ call flush(out) end if ! sending particles to down if (np_send_d.gt.0) then call MPI_SEND(xyzp1(1,np1-np_send_d+1),3*np_send_d,MPI_REAL8 ,id_down,6*myid+3,MPI_COMM_TASK,mpi_err) call MPI_SEND(uvwp1(1,np1-np_send_d+1),3*np_send_d,MPI_REAL8 ,id_down,6*myid+4,MPI_COMM_TASK,mpi_err) call MPI_SEND( ipart1(np1-np_send_d+1),np_send_d ,MPI_INTEGER4,id_down,6*myid+5,MPI_COMM_TASK,mpi_err) !!$ write(out,*) 'sent',np_send_d,' to the left:', ipart(np1-np_send_d+1:np1) !!$ call flush(out) end if ! completing the receive from up if (np_get_u.gt.0) then call MPI_WAIT(request1,mpi_status,mpi_err) call MPI_WAIT(request2,mpi_status,mpi_err) call MPI_WAIT(request3,mpi_status,mpi_err) !!$ write(out,*) 'received',np_get_u,' from the right:', ipart(np+1:np+np_get_u) !!$ call flush(out) np = np + np_get_u end if ! writing out particles' locations to separate files for each particle !!$ if (mod(ITIME,IWRITE4).eq.0) then if (nptot.lt.1000) then do i=1,np write(fname,"('p.',i4.4)") ipart(i) open(99,file=fname,position='append') write(99,'(i6.6,x, i4, 12e16.8)') & itime,myid,time,(xyzp(1:2,i)-1.0d0)*dx,(xyzp(3,i)-1.0d0+dble(myid*nz))*dz,& uvwp(:,i) close(99) end do end if !!$ end if return end subroutine particles_update_slabs !================================================================================ !================================================================================ !================================================================================ !================================================================================ !================================================================================ function cint_3d(a,x,y,z) result (blah) implicit none real*8 :: blah real*8 :: a(4,4,4), x, y, z logical :: compute_r integer :: j, k real*8 :: t2(4,4), t1(4), r(4) blah = 0.0d0 compute_r = .true. do k = 1,4 do j = 1,4 t2(j,k) = cint(a(1,j,k),a(2,j,k),a(3,j,k),a(4,j,k),x,r,compute_r) if (compute_r) compute_r = .false. end do end do compute_r = .true. do k = 1,4 t1(k) = cint(t2(1,k),t2(2,k),t2(3,k),t2(4,k),y,r,compute_r) if (compute_r) compute_r = .false. end do blah = cint(t1(1),t1(2),t1(3),t1(4),z,r,.true.) return end function cint_3d !================================================================================ function cint(p1,p2,p3,p4,x,r,compute_r) result (blah) implicit none logical :: compute_r real*8 :: p1,p2,p3,p4,x,r(4) real*8 :: sctmp real*8 :: blah ! If flag=.true., compute vector R(x) ! If flag=.false., consider vector R(x) defined and do not compute it if (compute_r) then sctmp = (x - 1.d0) * (x - 2.d0) r(1) = -x * sctmp r(2) = 3.d0 * (x + 1.d0) * sctmp sctmp = (x + 1) * x r(3) = -3.d0 * sctmp * (x - 2.d0) r(4) = sctmp * (x - 1.d0) end if blah = (r(1)*p1 + r(2)*p2 + r(3)*p3 + r(4)*p4) / 6.0d0 return end function cint end module m_particles