!====================================================================== ! ! MPI version ! ------------------------------------------- ! Purpose : Solve Ut+Ux+Uy=F(x,y) with Leap-Frog ! Author : Jarmo Rantakokko ! Date : 990615 ! !====================================================================== program adveq implicit none include 'mpif.h' integer, parameter :: DP=kind(0.0D0) !-- Variables. integer :: Nx,Ny,Nt,i,j,k,nnx,nny,rest,nx1,nx2,ny1,ny2 real(kind=DP) :: dt,dx,dy,norm,T,x,y,v,ti,sum real(kind=DP)::ttime,ctime,wtime,wtt real(kind=DP),pointer,dimension(:,:)::uold,u,unew,temp; integer,parameter:: disk=10 character(len=*),parameter :: input='params.dat' real(kind=DP)::F,up,h !-- MPI integer::ierror,myid,nproc,proc_grid,ndim integer,dimension(2)::coords,dims logical::reorder logical,dimension(2)::periods integer,pointer,dimension(:)::reqarr1,reqarr2,reqarr3,reqarrtemp integer::nreq,dummy integer,dimension(MPI_STATUS_SIZE,8)::mpistatus interface subroutine bound(u,t,coords,dims,nx1,ny1,dx,dy) integer, parameter :: DP=kind(0.0D0) real(kind=DP),pointer,dimension(:,:)::u real(kind=DP)::t,dx,dy integer::nx1,ny1 integer,dimension(2)::coords,dims end subroutine bound end interface !---------------------------------------------------------------------- ! Set up MPI call MPI_Init(ierror) reorder=.TRUE.;ndim=2;periods(1)=.FALSE.;periods(2)=.FALSE.; dims=0 call MPI_Comm_Size(MPI_COMM_WORLD, nproc, ierror) call MPI_Dims_create(nproc,ndim,dims,ierror) call MPI_Cart_create(MPI_COMM_WORLD,ndim,dims,periods,reorder, & proc_grid,ierror) call MPI_Comm_Rank(proc_grid, myid, ierror) call MPI_Cart_coords(proc_grid,myid,ndim,coords,ierror) ! Set up the problem namelist /problemsize/ Nx,Ny open(unit=disk,file=input) read(disk,problemsize) close(disk) ! Partition the data dx=1.0_DP/Nx; dy=1.0_DP/Ny; dt=1.0_DP/(Nx+Ny); T=1.0; Nt=nint(T/dt); nnx=(Nx+1)/dims(1); rest=mod((Nx+1),dims(1)); if (rest>coords(1)) then nnx=nnx+1; nx1=nnx*coords(1); nx2=nx1+nnx-1; else nx1=rest*(nnx+1)+(coords(1)-rest)*nnx; nx2=nx1+nnx-1; end if nny=(Ny+1)/dims(2); rest=mod((Ny+1),dims(2)); if (rest>coords(2)) then nny=nny+1; ny1=nny*coords(2); ny2=ny1+nny-1; else ny1=rest*(nny+1)+(coords(2)-rest)*nny; ny2=ny1+nny-1; end if ! Add ghost points (one on each side) allocate(uold(-1:nnx,-1:nny),u(-1:nnx,-1:nny)) allocate(unew(-1:nnx,-1:nny)); uold=0.0;u=0.0;unew=0.0 ! Initiate Persistent communication objects allocate(reqarr1(8),reqarr2(8),reqarr3(8)) call init_comm(uold,reqarr1,nreq,proc_grid,coords,dims) call init_comm(u,reqarr2,nreq,proc_grid,coords,dims) call init_comm(unew,reqarr3,nreq,proc_grid,coords,dims) if (myid.eq.0) then write(*,*) '============================================' write(*,*) 'Version : MPI-SunFire 15K' write(*,'(A,I8)') ' Number of nodes :',nproc write(*,'(A,3I8)') ' Problem Size :',Nx+1,Ny+1,Nt write(*,*) '============================================' write(*,*) 'Computing...' end if ! Total solve time call MPI_BARRIER(MPI_COMM_WORLD,ierror) ttime=MPI_WTIME() ! Initial conditions do j=0,nny-1 do i=0,nnx-1 x=real(i+nx1,kind=DP)/Nx; y=real(j+ny1,kind=DP)/Ny u(i,j)=h(x+y)+up(x,y); unew(i,j)=h(x+y-2*dt)+up(x,y); end do end do ! Integration time call MPI_BARRIER(MPI_COMM_WORLD,ierror) ctime=MPI_WTIME() wtime=0.0 ! Integrate the solution in time do k=2,Nt ! Swap pointers temp=>uold; uold=>u; u=>unew; unew=>temp; reqarrtemp=>reqarr1; reqarr1=>reqarr2 reqarr2=>reqarr3; reqarr3=>reqarrtemp ! Exchange partition boundary wtt=MPI_WTIME() call MPI_Startall(nreq,reqarr2,ierror) call MPI_Waitall(nreq,reqarr2,mpistatus,ierror) wtime=wtime+MPI_WTIME()-wtt ! Leap Frog do j=0,nny-1 do i=0,nnx-1 x=real(i+nx1,kind=DP)/Nx; y=real(j+ny1,kind=DP)/Ny unew(i,j)=uold(i,j)+2*dt*(F(x,y)- & ((u(i+1,j)-u(i-1,j))/2.0_DP*Nx+ & (u(i,j+1)-u(i,j-1))/2.0_DP*Ny)) end do end do ! Boundary conditions call bound(unew,k*dt,coords,dims,nx1,ny1,dx,dy) end do ! Integration time call MPI_BARRIER(MPI_COMM_WORLD,ierror) ctime=MPI_WTIME()-ctime ! Residual norm || u_new-(h(x+y-2*t)+up(x,y)) || sum=0.0 do j=0,nny-1 do i=0,nnx-1 x=real(i+nx1,kind=DP)/Nx; y=real(j+ny1,kind=DP)/Ny v=h(x+y-2*Nt*dt)+up(x,y) sum=sum+(unew(i,j)-v)*(unew(i,j)-v); end do end do norm=0.0 call MPI_Reduce(sum,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM,0, & proc_grid,ierror) ! Total solve time call MPI_BARRIER(MPI_COMM_WORLD,ierror) ttime=MPI_WTIME()-ttime ! Display results ! write(*,'(I3,F9.4,F9.4)') myid,ctime,wtime call MPI_BARRIER(MPI_COMM_WORLD,ierror) if (myid.eq.0) then write(*,*) '--------------' write(*,'(A,F10.4,A)') ' Solve time : ',ttime,' sec' write(*,'(A,E14.6)') ' Error norm : ',norm/sqrt(real(Nx*Ny)) write(*,*) '--------------' end if call MPI_Finalize(ierror) contains subroutine init_comm(u,reqarr,nreq,proc_grid,coords,dims) integer, parameter :: DP=kind(0.0D0) real(kind=DP),pointer,dimension(:,:)::u integer,dimension(:)::reqarr,coords,dims integer::nreq,proc_grid integer::left,right,up,down,ierror,status,STRIDED,Nx,Ny integer::tagx1=1,tagx2=2,tagy1=3,tagy2=4 Nx=ubound(u,1); Ny=ubound(u,2) call MPI_Cart_shift(proc_grid,0,1,left,right,ierror) call MPI_Cart_shift(proc_grid,1,1,up,down,ierror) call MPI_Type_vector(Ny,1,Nx+2,MPI_DOUBLE_PRECISION,STRIDED,ierror) call MPI_Type_commit(STRIDED,ierror) nreq=0 if (coords(1)>0) then nreq=nreq+1 call MPI_Send_init(u(0,0),1,STRIDED,left,tagx1,proc_grid, & reqarr(nreq),ierror) end if if (coords(1)0) then nreq=nreq+1 call MPI_Recv_init(u(-1,0),1,STRIDED,left,tagx2,proc_grid, & reqarr(nreq),ierror) end if if (coords(2)>0) then nreq=nreq+1 call MPI_Send_init(u(0,0),Nx+1,MPI_DOUBLE_PRECISION,up,tagy1, & proc_grid,reqarr(nreq),ierror) end if if (coords(2)0) then nreq=nreq+1 call MPI_Recv_init(u(0,-1),Nx+1,MPI_DOUBLE_PRECISION,up,tagy2, & proc_grid,reqarr(nreq),ierror) end if end subroutine init_comm end program adveq