PROGRAM MPIdif2Dcn ! An application to show MPI techniques applied to a mathematical problem ! Solve the 2D diffusion equation ! q_t = q_xx + q_yy ! on (x,y) in [0,1]x[0,1] with boundary conditions ! q(x,y,0)=sin(pi*x)*sin(pi*y) ! q(0,y,t)=q(1,y,t)=q(x,0,t)=q(x,1,t)=0 ! using the Crank-Nicolson algorithm. ! The computation domain is split into vertical strips. Each strip is ! sent to another process to carry out the time stepping. ! ! We solve the implicit system in the Crank-Nicolson algorithm using ! nGS Gauss-Seidel iterations per time step and a FTCS algorithm ! for the initial approximation to q at the new time step ! ! The linear system solved at each time step is ! ! (I - k/2 delta_h^2) qnew_ij = (I + k/2 delta_h^2) qold_ij = b_ij ! ! with delta_h^2 the discrete Laplace operator. ! ! The sequence of operations is ! ! (1) compute b from qold ! (2) use FTCS to get initial guess for qnew ! (3) carry out nGS Gauss-Seidel iterations ! (4) qold=qnew ! ! In this version we synchronize during the Gauss-Seidel iterations ! ! Compile with: ! mpif90 -o mpidif2Dcn mpidif2Dcn.f90 -lmpi ! Program requires data file mpidif2Dcn.dat. Sample data file: ! ! 64 64 ! 0.01 0.001 0.0001 ! 3 1.0 ! ! Run with: ! mpirun -c 2 -v mpidif2Dcn ! on a 2 processor computer or ! mpirun -v mpidif2Dcn ! on a cluster. ! Code prepared for Numerical Analysis Research Club (NARC) ! University of Washington ! Sorin Mitran, April, 2000 INCLUDE 'mpif.h' INTEGER, PARAMETER :: TAG_WRITEDATA=0, TAG_DATAWRITTEN=1 INTEGER, PARAMETER :: TAG_DATASYNCH=2, TAG_INIT=3 REAL, ALLOCATABLE, DIMENSION (:,:) :: qnew,qold REAL, ALLOCATABLE, DIMENSION (:) :: qleft,qright REAL, ALLOCATABLE, DIMENSION (:,:) :: b REAL qparam(7) DOUBLE PRECISION tIO,tcomp,tstart,tend,tplot1,tplot0 ! Initialize the MPI system CALL MPI_INIT(ierr) IF (ierr /= MPI_SUCCESS) THEN PRINT *,'Unable to initialize MPI system' STOP END IF ! Find out how many processes are active CALL MPI_COMM_SIZE(MPI_COMM_WORLD,np,ierr) ! Find out what the ID of this particular process is. CALL MPI_COMM_RANK(MPI_COMM_WORLD,idproc,ierr) IF (idproc == 0) THEN ! The master process reads input data OPEN(UNIT=1,FILE='mpidif2Dcn.dat',STATUS='OLD') READ(1,*)nx,ny READ(1,*)tfinal,tplot,delt READ(1,*)nGS,betaGS CLOSE(1) ! And then sends it to all the others qparam(1)=nx; qparam(2)=ny; qparam(3)=tfinal; qparam(4)=tplot; qparam(5)=delt qparam(6)=nGS; qparam(7)=betaGS DO n=1,np-1 CALL MPI_SEND(qparam(1), 7, MPI_REAL, n, TAG_INIT, MPI_COMM_WORLD, ierr) END DO ELSE ! Slave process receives input data from master CALL MPI_RECV(qparam(1), 7, MPI_REAL, 0, TAG_INIT, MPI_COMM_WORLD, ierr) nx=qparam(1); ny=qparam(2); tfinal=qparam(3); tplot=qparam(4); delt=qparam(5) nGS=qparam(6); betaGS=qparam(7) END IF ! Initialize. Each process allocates enough memory to hold the ! strip of data it's going to work on. ! The relationship between the local index i within a strip and the ! global index iglbl is ! iglbl=idproc*ni+i ! idproc is the process id running from 0 to np-1 ni=CEILING((nx+1.0)/np) ALLOCATE(qold(0:ni+1,0:ny),qnew(0:ni+1,0:ny),b(0:ni+1,0:ny),STAT=iError) ! Gridlines from 1 to ni are within this processor's strip. ! Gridline i=0 will be received from neighbor to the left ! (lower process number) or left boundary condition ! Gridline i=ni+1 will be received from neighbor to the right ! (higher process number) or right boundary condition IF (iError /= 0) THEN PRINT *,'Error allocating qold, qnew or b' STOP END IF ALLOCATE(qleft(0:ny),qright(0:ny),STAT=iError) IF (iError /= 0) THEN PRINT *,'Error allocating qleft or qright' STOP END IF ! Initial conditions. Since these do not depend on synchronization between ! the processes we can use the same code for all processes hx=1.0/nx; hy=1.0/ny; pi=4.0*atan(1.0); pi2=pi**2 DO j=0,ny piy=pi*j*hy DO i=0,ni+1 iglbl=idproc*ni+i pix=pi*iglbl*hx qold(i,j)=sin(pix)*sin(piy) qnew(i,j)=qold(i,j) END DO END DO ! Compute the time step to spatial step ratios sigmax=delt/hx**2; sigmay=delt/hy**2 sigmax2=0.5*sigmax; sigmay2=0.5*sigmay diagGS=1.0/(1.0+sigmax+sigmay) t=0.0; plottime=0.0 ! Start the time iteration. IF (idproc == 0) THEN PRINT *,'Running on ',np,' processes from t=0 to t=',tfinal,' with delt=',delt tcomp=0.0; tIO=0.0 END IF DO WHILE (t=plottime) THEN ! Time to output the current field plottime=plottime+tplot ! Set the next plot time ! There is no guaranteed running order among the processes. ! So we must set up a master/slave relationship in order to output ! the data in proper order. ! We also compute the norm-1 relative error at this time IF (idproc == 0) THEN qnorm=0.0; abserr=0.0 tplot0=MPI_WTIME() ! Master process. Open the output file PRINT *,'Plot time t=',t PRINT 2000,0 2000 FORMAT('Output of data from process ',i3) IF (t == 0.0) THEN OPEN(UNIT=2,FILE='mpidif2Dcn.out',STATUS='REPLACE') ELSE OPEN(UNIT=2,FILE='mpidif2Dcn.out',STATUS='OLD',POSITION='APPEND') END IF WRITE(2,2001)nx+1,ny+1 2001 FORMAT('ZONE I=',i4,' J=',i4) DO i=1,MIN(ni,nx+1-idproc*ni) iglbl=idproc*ni+i pix=pi*iglbl*hx DO j=0,ny piy=pi*j*hy WRITE(2,*)qnew(i,j) qexact=EXP(-2*pi2*t)*SIN(pix)*SIN(piy) abserr=abserr+ABS(qexact-qnew(i,j)) qnorm=qnorm+ABS(qexact) END DO END DO ! Close the output file. CLOSE(2) ! Loop over the processes telling each one to output its data DO n=1,np-1 PRINT 2000,n ! Send a request to output the data. CALL MPI_SEND(qsend, 1, MPI_REAL, n, TAG_WRITEDATA, MPI_COMM_WORLD, ierr) ! Wait for a response that data has been written. ! Also, receive error information. CALL MPI_RECV(qparam(1), 2, MPI_REAL, n, TAG_DATAWRITTEN, MPI_COMM_WORLD, ierr) abserr=abserr+qparam(1); qnorm=qnorm+qparam(2) END DO PRINT 1999,abserr/qnorm 1999 FORMAT('Relative error at this time eps=',e8.2) tplot1=MPI_WTIME() tIO=tIO+(tplot1-tplot0) ELSE ! Wait for a signal from the master process CALL MPI_RECV(qrecv, 1, MPI_REAL, 0, TAG_WRITEDATA, MPI_COMM_WORLD, ierr) abserr=0.0; qnorm=0.0 ! After signal is received, output the data OPEN(UNIT=2,FILE='mpidif2Dcn.out',STATUS='OLD',POSITION='APPEND') DO i=1,MIN(ni,nx+1-idproc*ni) iglbl=idproc*ni+i pix=pi*iglbl*hx DO j=0,ny piy=pi*j*hy WRITE(2,*)qnew(i,j) qexact=EXP(-2*pi2*t)*SIN(pix)*SIN(piy) abserr=abserr+ABS(qexact-qnew(i,j)) qnorm=qnorm+ABS(qexact) END DO END DO CLOSE(2) qparam(1)=abserr; qparam(2)=qnorm CALL MPI_SEND(qparam(1), 2, MPI_REAL, 0, TAG_DATAWRITTEN, MPI_COMM_WORLD, ierr) END IF END IF IF (idproc == 0) tstart=MPI_WTIME() ! Carry out a time step. All BC's are initially valid ! (1) compute rhs term b from qold DO j=1,ny-1 DO i=1,MIN(ni,nx+1-idproc*ni) b(i,j)=(1-sigmax-sigmay)*qold(i,j) & +sigmax2*(qold(i+1,j)+qold(i-1,j)) & +sigmay2*(qold(i,j+1)+qold(i,j-1)) END DO END DO ! (2) Use FTCS to compute initial approximation to qnew DO j=1,ny-1 DO i=1,MIN(ni,nx+1-idproc*ni) qnew(i,j)=qold(i,j)+sigmax*(qold(i+1,j)-2*qold(i,j)+qold(i-1,j)) & +sigmay*(qold(i,j+1)-2*qold(i,j)+qold(i,j-1)) END DO END DO ! (3) Carry out nGS Gauss-Seidel iterations on qnew DO n=1,nGS DO j=1,ny-1 DO i=1,MIN(ni,nx+1-idproc*ni) qnew(i,j)=diagGS*(b(i,j)+sigmax2*(qnew(i+1,j)+qnew(i-1,j)) & +sigmay2*(qnew(i,j+1)+qnew(i,j-1))) END DO END DO ! Synchronize across the processes. ! Send right ghost cell values IF (idproc < np-1) THEN qright(:)=qnew(ni,:) CALL MPI_SEND(qright(0), ny+1, MPI_REAL, idproc+1, TAG_DATASYNCH, MPI_COMM_WORLD, ierr) END IF ! Receive left ghost cell values IF (idproc > 0) THEN CALL MPI_RECV(qleft(0), ny+1, MPI_REAL, idproc-1, TAG_DATASYNCH, MPI_COMM_WORLD, ierr) qnew(0,:)=qleft(:) END IF ! Send left ghost cell values IF (idproc > 0) THEN qleft(:)=qnew(1,:) CALL MPI_SEND(qleft(0), ny+1, MPI_REAL, idproc-1, TAG_DATASYNCH, MPI_COMM_WORLD, ierr) END IF ! Receive right ghost cell values IF (idproc < np-1) THEN CALL MPI_RECV(qright(0), ny+1, MPI_REAL, idproc+1, TAG_DATASYNCH, MPI_COMM_WORLD, ierr) qnew(ni+1,:)=qright(:) END IF END DO IF (idproc == 0) THEN tend=MPI_WTIME() tcomp=tcomp+tend-tstart END IF ! Prepare for next iteration qold=qnew t=t+delt END DO IF (idproc == 0) THEN PRINT *,'Time spent in I/O t=',tIO PRINT *,'Computation time t=',tcomp END IF ! Shut down the MPI system CALL MPI_FINALIZE(ierr) END