!!==================================================================================== module particles !!==================================================================================== logical :: user_defines_particles=.TRUE. character(len=80) :: particle_step_flag real(kind=8),allocatable :: positions(:,:) real(kind=8),allocatable :: uvels (:,:) real(kind=8),allocatable :: vvels (:,:) real(kind=8),allocatable :: wvels (:,:) integer :: nparticles integer :: j_particle_time integer :: istart_trajectories=0 !!==================================================================================== end module particles !!==================================================================================== subroutine InitializeParticles use user_routines use particles use user_data use mpi_params use etc, only: docfile implicit none integer :: i, j, npts real(kind=8) :: urv real(kind=8) :: sides(3) call random_seed !! I'm using the f90 intrinsic functions !! random_seed, random_number here j_particle_time = 1 npts = nparticles allocate( positions(npts,3) ) allocate( uvels(npts,4) ) !! 4 time levels allocate( vvels(npts,4) ) allocate( wvels(npts,4) ) sides(:)=(/Lx,Ly,Lz/) user_defines_particles = .TRUE. if (user_defines_particles) then call user_particle_positions(positions,npts,Lx,Ly,Lz) else !! I'll distribute them uniformly do i=1,npts do j=1,3 call RANDOM_NUMBER( urv ) positions(i,j) = sides(j)*urv enddo enddo endif positions=positions/Lz uvels(:,:) = 0.d0 vvels(:,:) = 0.d0 wvels(:,:) = 0.d0 if(myid==0) then open(1,file=docfile,position='append') write(1,*) ' -----> InitializeParticles routine exiting normally <---------- ' write(0,*) ' -----> InitializeParticles routine exiting normally <---------- ' close(1) endif end subroutine InitializeParticles !============================================================================ !============================================================================ !=============================================================== !=============================================================== subroutine advance_particles use etc, only: PM0,PM1,PM2,PM3, & my_z_limits,t_secs, & numzplanes_below,M_oldest use dependent_variables, only: velocity use intermediate_variables, only: explicit_rhs use user_data, only: dt,U0,Lx,Ly,Lz, & xdim_periodic, & ydim_periodic, & zdim_periodic,pi,ny use independent_variables use particles use mpi_params implicit none include 'mpif.h' integer :: i,j,k,kk,ngrdpts,ierr real(kind=8) :: xx,yy,zz,ans1,ans2,ans3 integer, SAVE :: m(3) integer, allocatable,SAVE :: my_labels(:) real(kind=8),allocatable,SAVE :: tmp(:) real(kind=8), save :: dt2,dt12,dt24 logical, save :: first_entry = .TRUE. if(nparticles<=0) return !!================================================================= !!================================================================= if( first_entry ) then dt2=dt/2.d0 dt12=dt/12.d0 dt24=dt/24.d0 m(:)=0 !! interpolate function vals, not derivs allocate( my_labels(nparticles),tmp(nparticles) ) first_entry = .FALSE. endif !!================================================================= !!================================================================= !!================================================================= !!================================================================= !! periodic wrapping etc !!================================================================= !!================================================================= if(xdim_periodic) then do i=1,nparticles if( positions(i,1) < 0 ) then positions(i,1) = positions(i,1) + Lx/Lz elseif( positions(i,1) > Lx/Lz ) then positions(i,1) = positions(i,1) - Lx/Lz endif enddo endif if(ydim_periodic .and. ny>1) then do i=1,nparticles if( positions(i,2) < 0 ) then positions(i,2) = positions(i,2) + Ly/Lz elseif( positions(i,2) > Ly/Lz ) then positions(i,2) = positions(i,2) - Ly/Lz endif enddo endif if(zdim_periodic) then do i=1,nparticles if( positions(i,3) < 0 ) then positions(i,3) = positions(i,3) + Lz/Lz elseif( positions(i,3) > Lz/Lz ) then positions(i,3) = positions(i,3) - Lz/Lz endif enddo endif !!================================================================= !!================================================================= !!================================================================= !! mark the indices of particles residing on processor myid !!================================================================= do i=1,nparticles if( positions(i,3) >= my_z_limits(1) .and. & positions(i,3) < my_z_limits(2) ) then my_labels(i) = 1 else my_labels(i) = 0 positions(i,1:3) = 0.d0 endif enddo !!================================================================= !!================================================================= !!================================================================= !!================================================================= !! get the u,v,w values at the positions of the particles on myid !! arrays are dimensioned like !! positions(npts,3) !! 3 space dims !! uvels(npts,4) !! 4 time levels !!================================================================= !!================================================================= !! get interpolated u values ( explicit_rhs are scratch arrays ) j=M_oldest !! oldest of the MM0,MM1,MM2,MM3 set, can be trashed call spline_interp_value(positions(1,1), & positions(1,2), & positions(1,3), & nparticles, & velocity(1), & m(1), & uvels(1,PM0), & my_labels, & explicit_rhs(1,j)%SlabVals, & explicit_rhs(2,j)%SlabVals, & explicit_rhs(3,j)%SlabVals ) !! get interpolated v values call spline_interp_value(positions(1,1), & positions(1,2), & positions(1,3), & nparticles, & velocity(2), & m(1), & vvels(1,PM0), & my_labels, & explicit_rhs(1,j)%SlabVals, & explicit_rhs(2,j)%SlabVals, & explicit_rhs(3,j)%SlabVals ) !! get interpolated w values call spline_interp_value(positions(1,1), & positions(1,2), & positions(1,3), & nparticles, & velocity(3), & m(1), & wvels(1,PM0), & my_labels, & explicit_rhs(1,j)%SlabVals, & explicit_rhs(2,j)%SlabVals, & explicit_rhs(3,j)%SlabVals ) !!================================================================= !!================================================================= !! At this point positions(i,:) are nonzero only if !! particle i resides on processor myid. Similarly !! for uvels(i,PM0),vvels(i,PM0) etc. Values for !! PM1, PM2 etc should be identical for all processors. !!================================================================= !!================================================================= !! Integrate forward one time step. do i=1,nparticles if( my_labels(i)==1 ) then !! only deal with those on proc. myid if( particle_step_flag .eq. 'euler' ) then positions(i,1)=positions(i,1) + dt*uvels(i,PM0) positions(i,2)=positions(i,2) + dt*vvels(i,PM0) positions(i,3)=positions(i,3) + dt*wvels(i,PM0) particle_step_flag='AB2' elseif( particle_step_flag .eq. 'AB2' ) then positions(i,1)=positions(i,1) + dt2*( 3.*uvels(i,PM0) & - uvels(i,PM1) ) positions(i,2)=positions(i,2) + dt2*( 3.*vvels(i,PM0) & - vvels(i,PM1) ) positions(i,3)=positions(i,3) + dt2*( 3.*wvels(i,PM0) & - wvels(i,PM1) ) particle_step_flag='AB3' elseif( particle_step_flag .eq. 'AB3' ) then positions(i,1)=positions(i,1) + dt12*(23.*uvels(i,PM0) & - 16.*uvels(i,PM1) & + 5.*uvels(i,PM2) ) positions(i,2)=positions(i,2) + dt12*(23.*vvels(i,PM0) & - 16.*vvels(i,PM1) & + 5.*vvels(i,PM2) ) positions(i,3)=positions(i,3) + dt12*(23.*wvels(i,PM0) & - 16.*wvels(i,PM1) & + 5.*wvels(i,PM2) ) particle_step_flag='AB4' elseif( particle_step_flag .eq. 'AB4' ) then positions(i,1)=positions(i,1) + dt24*(55.*uvels(i,PM0) & - 59.*uvels(i,PM1) & + 37.*uvels(i,PM2) & - 9.*uvels(i,PM3) ) positions(i,2)=positions(i,2) + dt24*(55.*vvels(i,PM0) & - 59.*vvels(i,PM1) & + 37.*vvels(i,PM2) & - 9.*vvels(i,PM3) ) positions(i,3)=positions(i,3) + dt24*(55.*wvels(i,PM0) & - 59.*wvels(i,PM1) & + 37.*wvels(i,PM2) & - 9.*wvels(i,PM3) ) endif endif enddo !!================================================================= !!================================================================= !! Now globally_update_particle_data !!================================================================= !!================================================================= tmp(:) = uvels(:,PM0) call mpi_allreduce(tmp,uvels(1,PM0),nparticles, & mpi_double_precision,mpi_sum,comm,ierr ) tmp(:) = vvels(:,PM0) call mpi_allreduce(tmp,vvels(1,PM0),nparticles, & mpi_double_precision,mpi_sum,comm,ierr ) tmp(:) = wvels(:,PM0) call mpi_allreduce(tmp,wvels(1,PM0),nparticles, & mpi_double_precision,mpi_sum,comm,ierr ) do i=1,3 tmp(:) = positions(:,i) call mpi_allreduce(tmp,positions(1,i),nparticles, & mpi_double_precision,mpi_sum,comm,ierr ) enddo !!================================================================= !!================================================================= !! All processors now have updated positions and velocities !! for all particles, whether they reside locally or not !!================================================================= !!================================================================= end subroutine advance_particles !=============================================================== !=============================================================== subroutine spline_interp_value(xvec,yvec,zvec,npts,FF, & m,ans,my_labels, & s_x,s_y,s_z, & nhat_1,nhat_2,nhat_3) !! !! (1) Given a vector of (dimensionless) x-y-z locations !! calculate the corresponding parameter values s. !! e.g. x=f(s) given x* find s* such that f(s*)=x* !! Do this in each coordinate direction x-y-z. !! !! (2) Call "spline" to get the interpolated field value !! at each point. This part of the problem is done in !! the space {s1, s2, s3}, the equally spaced parameters !! for coordinates x,y,z. The s* values found in (1) define !! the current particle positions in this space. !! if !! m(i) = 0,1,2 return F, F' or F'' wrt the "ith" dimension !!============================================================= use class_ScalarVariable use user_routines use user_data, only: Lx,Ly,Lz, & nx,ny,nz, & xdim_periodic, & ydim_periodic, & zdim_periodic implicit none type( ScalarVariable ) :: FF integer :: m(3) integer :: npts,i integer :: my_labels(npts) real(kind=8) :: xvec(npts) real(kind=8) :: yvec(npts) real(kind=8) :: zvec(npts) real(kind=8) :: s_x (npts) real(kind=8) :: s_y (npts) real(kind=8) :: s_z (npts) real(kind=8),optional :: nhat_1(npts) real(kind=8),optional :: nhat_2(npts) real(kind=8),optional :: nhat_3(npts) real(kind=8) :: locs(3) real(kind=8) :: svals(3) real(kind=8) :: ans(npts) !! Get the parametric (s) values corresponding to xyz positions. do i=1,npts if( my_labels(i)==1 ) then !! get dimensional x,y,z for trajectory i locs(:)=Lz*(/xvec(i),yvec(i),zvec(i)/) !! get the corresponding parameter values s_x,s_y,s_z call point_map_x2s(locs(1),1,xdim_periodic,Lx,nx,s_x(i)) call point_map_x2s(locs(2),2,ydim_periodic,Ly,ny,s_y(i)) call point_map_x2s(locs(3),3,zdim_periodic,Lz,nz,s_z(i)) endif enddo if( present(nhat_1) .and. present(nhat_2) .and. present(nhat_3) ) then call spline(FF,ans,s_x,s_y,s_z,m,npts,my_labels,nhat_1,nhat_2,nhat_3) else call spline(FF,ans,s_x,s_y,s_z,m,npts,my_labels) endif end subroutine spline_interp_value !=============================================================== !=============================================================== !!=============================================================== !!=============================================================== subroutine spline(F,g,x,y,z,m,npts,my_labels,nhat_1,nhat_2,nhat_3) !! This function interpolates three-dimensional discrete !! data to a set of arbitary, i.e. non-gridpoint, !! locations x(1:npts),y(1:npts),z(1:npts). !! given the discrete gridpoint values F(:,:,:). !! Results are returned in g(1:npts). !! !! m(idim) return m-th derivative, in the idim direction, of F. !! m can be 0, 1 or 2. !! !! Required subroutines !! db3ink: computes coefficients of the 3D tensor basis functions !! db3val: evaluates the tensor product spline !! db2ink: computes coefficients of the 2D tensor basis functions !! db2val: evaluates the correspondingtensor product spline !! !! routine requires no interprocessor communication !! but, rather, uses the 'not-a-knot' end conditions !! and computes near edge values locally. Similarly, !! the spline basis does not explicitly impose periodicity, !! again relying on the 'not-a-knot' end condition. !!=============================================================== !!=============================================================== use Class_ScalarVariable use mpi_params use independent_variables use user_routines, only: point_map_s2x use user_data, only: Lx,Ly,Lz, & xdim_periodic, & ydim_periodic, & zdim_periodic use etc, only: numzplanes, & numzplanes_below implicit none type( ScalarVariable ) :: F integer :: i,j,k integer,SAVE :: n1,n2,n3 integer,SAVE :: nx,ny,nz,globalnz integer :: m(3),npts,k0,k1 integer,SAVE :: p,order(3),nwrk,iflag,idim integer :: my_labels(npts) real(kind=8) :: x(npts),xval real(kind=8) :: y(npts),yval real(kind=8) :: z(npts),zval real(kind=8) :: g(npts) real(kind=8),allocatable,SAVE :: Fext(:,:,:) real(kind=8),allocatable,SAVE :: xgrid(:) real(kind=8),allocatable,SAVE :: ygrid(:) real(kind=8),allocatable,SAVE :: zgrid(:) real(kind=8),allocatable,SAVE :: tx(:) real(kind=8),allocatable,SAVE :: ty(:) real(kind=8),allocatable,SAVE :: tz(:) real(kind=8),allocatable,SAVE :: e3(:,:,:) real(kind=8),allocatable,SAVE :: wrk(:) real(kind=8),optional :: nhat_1(npts) real(kind=8),optional :: nhat_2(npts) real(kind=8),optional :: nhat_3(npts) real(kind=8),external :: db3val,db2val real(kind=8) :: dx,dy,dz,h real(kind=8) :: f_s,f_x,f_y,f_z,x_s,y_s,z_s real(kind=8) :: Fmin,Fmax logical,SAVE :: periodic(3) logical,SAVE :: limits_test logical,SAVE :: first_entry=.TRUE. if( first_entry ) then limits_test = .FALSE. p=3 ! cubic order(:)=p+1 periodic=(/xdim_periodic, & ydim_periodic, & zdim_periodic/) call GetLocalArraySize(F,n1,n2,n3) !! i.e. nx,ny,locnz for F !! determine the size of the data array we'll actually work with if(xdim_periodic) then nx=n1+1 else nx=n1 endif if(ydim_periodic .and. n2>1) then ny=n2+1 else ny=n2 endif if(zdim_periodic .and. myid==numprocs-1) then nz=n3+1 elseif( .not. zdim_periodic .and. myid==numprocs-1 ) then nz=n3 endif if( myid1) then h = 1./(ny-1.) ygrid(1:ny)=(/( (i-1.)*h,i=1,ny )/) elseif(ny==1) then ygrid(1)=1.d0 endif !! these are the s positions on myid k0 = numzplanes_below(myid) + 1 k1 = k0 + numzplanes(myid) - 1 zgrid(1:n3) = Grid(0)%SpaceDims(3)%s(k0:k1) globalnz = ubound(Grid(0)%SpaceDims(3)%s,1) !! if we are extending myid's data upward, get the next value as well if(nz>n3) zgrid(nz) = Grid(0)%SpaceDims(3)%s(k1+1) iflag=0 !! use default knot selection scheme if( order(1) .ge. nx ) order(1)=nx-1 if( order(2) .ge. ny ) order(2)=ny-1 if( order(3) .ge. nz ) order(3)=nz-1 first_entry = .FALSE. endif !! fill extended array call AssignFextValues(F,n1,n2,n3,Fext,nx,ny,nz,periodic) !!======================================================== !! Compute coefficients of the 3D tensor basis functions !! using deBoor routine from nist/gams/cmlib !!======================================================== if(ny>1) then call db3ink(xgrid,nx,ygrid,ny,zgrid,nz, & Fext(1,1,1),nx,ny, & order(1),order(2),order(3), & tx(1),ty(1),tz(1), & e3,wrk,iflag) if(iflag .ne. 1 ) stop 'error condition detected by db3ink' elseif(ny==1) then call db2ink(xgrid,nx,zgrid,nz,Fext,nx,order(1),order(3), & tx(1),tz(1),e3,wrk,iflag) if(iflag .ne. 1 ) stop 'error condition detected by db2ink' endif !!=========================================================== !! Now loop through each position & interpolate the function !!=========================================================== !!=========================================================== !! CASE 1: m(:)=0 evaluate function only, !! and only when my_labels = 1 !! do 2d ny=1 case separately !!=========================================================== if( m(1)==0 .and. m(2)==0 .and. m(3)==0 ) then if(ny>1) then do i=1,npts !! note, positions detected as outside the appropriate !! range will trigger a fatal error if( my_labels(i)==1 ) then g(i) = db3val(x(i),y(i),z(i), & m(1),m(2),m(3), & tx,ty,tz, & nx,ny,nz, & order(1),order(2),order(3), & e3,wrk) else g(i)=0.d0 endif enddo elseif(ny==1) then do i=1,npts !! note, positions detected as outside the appropriate !! range will trigger a fatal error if( my_labels(i)==1 ) then g(i) = db2val(x(i),z(i),m(1),m(3),tx,tz,nx,nz, & order(1),order(3),e3,wrk) else g(i)=0.d0 endif enddo endif endif !!=========================================================== !! CASE 2: m(1)=0 evaluate Fx,Fy,Fz and dot result with n_hat !! NB don't bother with my_labels here !! handle 2d case as above !! Also, note that values returned from point_map_s2x are !! DIMENSIONAL and so need to be scaled by Lz !! Also, though the inputs to point_map_s2x are stored in !! arrays called x,y,z they are really the s values in !! each of the 3 dimensions. !! !! {x_s}=x_s/Lz df/dx = {df/ds} * {ds/dx} = {df/ds}/{x_s} !! etc. !!=========================================================== if( m(1)==1 .and. m(2)==1 .and. m(3)==1 .and. present(nhat_1) & .and. present(nhat_2) .and. present(nhat_3) ) then if(ny>1) then do i=1,npts idim = 1 call point_map_s2x(x(i),idim,periodic(idim),Lx,n1,xval,x_s) f_s = db3val(x(i),y(i),z(i), & 1,0,0, & tx,ty,tz, & nx,ny,nz, & order(1),order(2),order(3), & e3,wrk) f_x = f_s / (x_s/Lz) idim = 2 call point_map_s2x(y(i),idim,periodic(idim),Ly,n2,yval,y_s) f_s = db3val(x(i),y(i),z(i), & 0,1,0, & tx,ty,tz, & nx,ny,nz, & order(1),order(2),order(3), & e3,wrk) f_y = f_s/(y_s/Lz) idim = 3 call point_map_s2x(z(i),idim,periodic(idim),Lz,globalnz,zval,z_s) f_s = db3val(x(i),y(i),z(i), & 0,0,1, & tx,ty,tz, & nx,ny,nz, & order(1),order(2),order(3), & e3,wrk) f_z = f_s/(z_s/Lz) !! dot product g(i) = nhat_1(i)*f_x + nhat_2(i)*f_y + nhat_3(i)*f_z enddo elseif(ny==1) then do i=1,npts idim = 1 call point_map_s2x(x(i),idim,periodic(idim),Lx,n1,xval,x_s) f_s = db2val(x(i),z(i),1,0,tx,tz,nx,nz, & order(1),order(3),e3,wrk) f_x = f_s/(x_s/Lz) idim = 3 call point_map_s2x(z(i),idim,periodic(idim),Lz,globalnz,zval,z_s) f_s = db2val(x(i),z(i),0,1,tx,tz,nx,nz, & order(1),order(3),e3,wrk) f_z = f_s/(z_s/Lz) !! dot product g(i) = nhat_1(i)*f_x + nhat_3(i)*f_z enddo endif endif end subroutine spline !=============================================================== !=============================================================== !=============================================================== !=============================================================== subroutine AssignFextValues(F,n1,n2,n3,Fext,nx,ny,nz,periodic) use class_ScalarVariable use mpi_params implicit none include 'mpif.h' type( ScalarVariable ) :: F logical :: periodic(3) integer :: n1,n2,n3 !! actual local array size of F%SlabVals integer :: nx,ny,nz !! size of extended array Fext integer :: i,j,k integer :: ierr,reqs(2) integer :: status_array(mpi_status_size,2) real(kind=8) :: Fext(nx,ny,nz) real(kind=8),allocatable,SAVE :: tmp(:,:) integer, SAVE :: src,count logical, SAVE :: first_entry=.TRUE. if( first_entry ) then allocate( tmp(n1,n2) ) first_entry = .FALSE. endif nreqs=0 status_array=0 src=myid+1 dest=myid-1 tag=99 count=n1*n2 if(myid /= numprocs-1) then !! all but top processor get ghost data from above ! recv 1 plane from above nreqs=nreqs+1 call MPI_IRECV(tmp,count,mpi_double_precision,src,tag,comm,reqs(nreqs),ierr) if( ierr .ne. 0 ) then stop 'MPI IRECV ERROR IN AssignFextValues/particle_routines.f90' endif endif !call mpi_barrier(comm,ierr) !! all recvs posted before sends begin ! Make these ISENDs, don't worry about buffering (2d) & eliminate the barrier if( myid /= 0 ) then !! all but bottom processor send ghost data below ! send bottom plane of F%SlabVals to myid-1 nreqs=nreqs+1 call MPI_ISEND(F%SlabVals,count,mpi_double_precision,dest,tag,comm,reqs(nreqs),ierr) if( ierr .ne. 0 ) then stop 'MPI IRSEND ERROR IN AssignFextValues/particle_routines.f90' endif endif !!================================================== !! do these copies while waiting for !! the messages to complete !!================================================== Fext(1:n1,1:n2,1:n3) = F%SlabVals(1:n1,1:n2,1:n3) if(periodic(1)) then Fext(nx,1:n2,1:n3)=F%SlabVals(1,1:n2,1:n3) endif if(periodic(2) .and. n2 > 1) then Fext(1:n1,ny,1:n3)=F%SlabVals(1:n1,1,1:n3) endif if(periodic(1) .and. periodic(2)) then Fext(nx,ny,1:n3)=F%SlabVals(1,1,1:n3) endif !!================================================== !!make sure all myid's communications are done before proceeding if(nreqs>0) then call mpi_waitall(nreqs,reqs(1),status_array(1,1),ierr) endif !!================================================== if(nz>n3) then Fext(1:n1,1:n2,nz)=tmp(1:n1,1:n2) Fext(nx,:,nz)=Fext(1,:,nz) Fext(:,ny,nz)=Fext(:,1,nz) Fext(nx,ny,nz)=Fext(1,1,nz) endif end subroutine AssignFextValues