program trMPI c Solve reaction-diffusion in 3D with 1D domain decomposition c Ct=D*(Cxx+Cyy+Czz)-M(C) c Include, declare variables... include 'mpif.h' implicit real*8 (a-h,o-z) parameter(nx=84,ny=80,nz=200) integer kbeg,kend,mbbeg,mbend,ncbeg,ncend real*8 c0(0:nx+1,0:ny+1,0:nz+1),c1(0:nx+1,0:ny+1,0:nz+1) parameter(mpi_master=0) c etc... c MPI Init call MPI_INIT( mpi_err ) call MPI_COMM_SIZE( MPI_COMM_WORLD, numnodes, mpi_err ) call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpi_err ) t1mpi=MPI_Wtime() c Divide domain into numnodes slabs in z-direction nzmpi=nz/numnodes dzmpi=float(nz)/float(numnodes) kbeg=int(float(myid)*dzmpi+1.e0) kend=int(1.e0+float(myid+1)*dzmpi)-1 kbeg2=kbeg if (kbeg.eq.1) kbeg2=2 kend2=kend if (kend.eq.nz) kend2=nz+1 kbeg3=kbeg if (kbeg.eq.1) kbeg3=0 myrow=myid/numnodes+1 mycol=1 c Make the row and col communicators c all processors with the same row will be in the same ROW_COMM call MPI_COMM_SPLIT(MPI_COMM_WORLD,myrow,mycol,ROW_COMM,mpi_err) call MPI_COMM_RANK( ROW_COMM, myid_row, mpi_err ) call MPI_COMM_SIZE( ROW_COMM, nodes_row, mpi_err ) c Find id of neighbors using the communicators created above myleft = myid_row-1;if( myleft .lt. 0 )myleft 1 = MPI_PROC_NULL myright = myid_row+1;if( myright .eq. nodes_row )myright 1 = MPI_PROC_NULL c Set parameters... c Open data, result files open(24,file= 'tran22G.su24a2Ai') open(10,file= 'cot33SPC.a2Aiy') cofile='conT33SPC.a2Aiy' cifile=cofile c MPI Variable Read/Init (master node) c Initialize variables or read in continuation data if(myid .eq. mpi_master) then open(15,file=cifile) if (ncont.eq.1) then read(15,*)ti0,tf0,nx0,ny0,nz0 else nx0=nx ny0=ny nz0=nz endif do 350 i0=1,nx0+1 do 350 j0=1,ny0+1 do 350 k0=1,nz0+1 if (ncont.eq.1) then read(15,*)c0(i,j,k) else c0(i,j,k)=1.d-5*(1.d0+dcos(xpi*dfloat(k-1)/dfloat(nz+1))) endif 350 continue close(15) endif c Broadcast data to all nodes call MPI_BCAST(c0(0:nx+1,0:ny+1,0:nz+1),(nx+2)*(ny+2)*(nz+2), 1 MPI_DOUBLE_PRECISION,mpi_master, 2 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(sseg0,nc,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) c MPI Geometry Read/Init (master node) if(myid .eq. mpi_master) then read(24,*)xl,yl,zl do 10 k=1,nc read(24,*) x1(k),y1(k),z1(k) read(24,*) ninfl(k),minfl(k,1),minfl(k,2) 10 continue close(24) endif c Broadcast geometry data call MPI_BCAST(xl,1,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(yl,1,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(zl,1,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(x1,nc,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(y1,nc,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(z1,nc,MPI_DOUBLE_PRECISION,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(ninfl,nc,MPI_INTEGER,mpi_master, 1 MPI_COMM_WORLD,mpi_err) call MPI_BCAST(minfl(1:nc,1:2),2*nc,MPI_INTEGER,mpi_master, 1 MPI_COMM_WORLD,mpi_err) c Time-stepping loop dt=1.0d-3 ti=0.0d0 tf=30.0d0 nt=int((tf-ti+1.d-8)/dt)+1 do 99 n=1,nt tc=ti+dt*dfloat(n-1) c Output data every nskip steps for point on 'end' node if (kend.eq.nz) then if(mod(n-1,nskip).eq.0) then write(10,*)sngl(tc),sngl(c0(nx,ny,nz)/sol) endif endif c Advance solution at internal points 22 do 120 k=kbeg2,kend do 120 j=istart,ny do 120 i=istart,nx mc=mc0v(i,j,k)*c0(i,j,k)/(ccrit+c0(i,j,k)) dctmp=dti*(rhx*(c0(i-1,j,k) 1 +c0(i+1,j,k)-2.d0*c0(i,j,k))+rhy*(c0(i,j-1,k) 2 +c0(i,j+1,k)-2.d0*c0(i,j,k))+rhz*(c0(i,j,k-1) 3 +c0(i,j,k+1)-2.d0*c0(i,j,k)))-dt*mc c1(i,j,k)=c0(i,j,k)+dctmp 120 continue c Apply BC: periodic in x and y, no flux in z do 156 k=kbeg2,kend do 156 i=1,nx c1(i,ny+1,k)=c1(i,1,k) c1(i,0,k)=c1(i,ny,k) 156 continue do 157 k=kbeg2,kend do 157 j=1,ny c1(nx+1,j,k)=c1(1,j,k) c1(0,j,k)=c1(nx,j,k) 157 continue do 158 k=kbeg2,kend c1(nx+1,ny+1,k)=c1(1,1,k) c1(nx+1,0,k)=c1(1,ny,k) c1(0,ny+1,k)=c1(nx,1,k) c1(0,0,k)=c1(nx,ny,k) 158 continue do 159 i=0,nx+1 do 159 j=0,ny+1 if (kend.eq.nz) then c1(i,j,nz+1)=0.d0 if(ind(i,j,nz).ne.-1) c1(i,j,nz+1)=c1(i,j,nz) if(ind(i,j,nz-1).ne.-1) then c1(i,j,nz+1)=(4.d0*c1(i,j,nz)-c1(i,j,nz-1))/3.d0 endif endif if (kbeg2.eq.2) then c1(i,j,1)=0.d0 if(ind(i,j,2).ne.-1) c1(i,j,1)=c1(i,j,2) if(ind(i,j,3).ne.-1) then c1(i,j,1)=(4.d0*c1(i,j,2)-c1(i,j,3))/3.d0 endif endif 159 continue c Update solution do 200 k=kbeg,kend2 do 200 j=istart-1,ny+1 do 200 i=istart-1,nx+1 c0(i,j,k)=c1(i,j,k) 200 continue c Exchange data on ends of slabs ireq=0 if(myleft .ne. MPI_PROC_NULL)then c send to left ireq=ireq+1 sv3=c0(:,:,kbeg) call MPI_ISEND(sv3,(nx+2)*(ny+2) ,MPI_DOUBLE_PRECISION,myleft, 1 100,ROW_COMM,req(ireq),mpi_err) c rec from left ireq=ireq+1 call MPI_IRECV(rv3,(nx+2)*(ny+2),MPI_DOUBLE_PRECISION,myleft, 1 100,ROW_COMM,req(ireq),mpi_err) endif if(myright .ne. MPI_PROC_NULL)then c rec from right ireq=ireq+1 call MPI_IRECV(rv4,(nx+2)*(ny+2),MPI_DOUBLE_PRECISION,myright, 1 100,ROW_COMM,req(ireq),mpi_err) c send to right ireq=ireq+1 sv4=c0(:,:,kend) call MPI_ISEND(sv4,(nx+2)*(ny+2),MPI_DOUBLE_PRECISION,myright, 1 100,ROW_COMM,req(ireq),mpi_err) endif call MPI_WAITALL(ireq,req,stat_ray,mpi_err) if (myleft .ne. MPI_PROC_NULL) c0(:,:,kbeg-1)=rv3 if (myright .ne. MPI_PROC_NULL) c0(:,:,kend+1)=rv4 99 continue c End time-stepping loop c Gather all results to master node for output do i=0,nx+1 do j=0,ny+1 call MPI_GATHERV(c0(i,j,kbeg3:kend2), 1 kend2-kbeg3+1,MPI_DOUBLE_PRECISION, 2 c02(i,j,:),counts, 4 offsets,MPI_DOUBLE_PRECISION, 3 mpi_master,MPI_COMM_WORLD,mpi_err) enddo enddo c Final output, close files if (myid.eq.mpi_master) then open(15,file=cofile) write(15,*)ti,tf,nx,ny,nz do 352 i=1,nx+1 do 352 j=1,ny+1 do 352 k=1,nz+1 write(15,*)c02(i,j,k) 352 continue endif close(15) close(10) endif c Calculate, output runtime t2mpi=MPI_Wtime() call MPI_REDUCE(t2mpi,t2mpit,1,MPI_DOUBLE_PRECISION, 1 MPI_SUM,mpi_master,MPI_COMM_WORLD,mpi_err) call MPI_REDUCE(t1mpi,t1mpit,1,MPI_DOUBLE_PRECISION, 1 MPI_SUM,mpi_master,MPI_COMM_WORLD,mpi_err) if (myid.eq.mpi_master) write(6,*)'runtime=',t2mpi-t1mpi if (myid.eq.mpi_master) write(6,*)'runtot=',t2mpit-t1mpit C End MPI call MPI_Finalize(mpi_err) stop end