subroutine etchn(iam, nodes, ierr) c- include 'mpif.h' c integer iam, nodes integer ierr, status(MPI_STATUS_SIZE) c double precision stime, etime, htime double precision sinit,endtime,wholtime real*8 r,dt,dx,dy,w,tol,tol1,t,b,r2,r3 real*8 d2x,d2y,d2x2,d2y2,deltay real*8 u1t(1650,1650),u1new(1650,1650) real*8 u1old(1650,1650),tmp(1650) real*8 u2t(825,1650),u2new(825,1650),u2old(825,1650), 1 temp1(1650,1650) c integer ctime,cstime,cetime,itercnt,time,ind1,ind2 integer itercnt,ind1,ind2, idummy1, idummy2 double precision ctime,cstime,cetime,time, dummy1, dummy2 integer maxrow1,maxcol1,maxrow2,maxcol2,a1,a2,i,j,l,maxsiz integer sum,flag,done,temp,maxtime,itime,row integer allnodes,npid,host,internod,maxnod integer intsiz,realsiz,matsiz,vecsiz,lastrow,transiz,transiz1 integer numrows1,numrows2,extra1,extra2 c- common /both/ dt,dx,dy,w,tol,b,r common /ireg1/ maxrow1,maxcol1,a1,maxrow2,maxcol2,a2 common /inidt1/ maxnod,internod,maxtime common /inidt2/ t common /top/temp1 common /bottom/u2t common /newtop/u1new common /newbot/u2new c common /perf/time,ctime,itercnt common /perf/time,ctime,dummy1 common /iperf/idummy1, idummy2, itercnt common /rows/numrows1,numrows2,extra1,extra2 c parameter( maxsiz = 1650) parameter( matsiz = 2722500) parameter( vecsiz = 1650) parameter( transiz = 2722500) parameter(transiz1 = 1361250) c c matsiz = maxsiz*maxsiz c vecsiz = maxsiz c if (iam .ne. 0) call MPI_RECV(maxnod ,3, MPI_INTEGER, 0, 1 51, MPI_COMM_WORLD, status, ierr ) if (iam .ne. 0) call MPI_RECV(t, 1, MPI_DOUBLE_PRECISION, 0, 1 55, MPI_COMM_WORLD, status, ierr ) if (iam .ne. 0) call MPI_RECV(numrows1, 4, MPI_INTEGER, 0, 1 75, MPI_COMM_WORLD, status, ierr ) if (iam .ne. 0) call MPI_RECV(dt, 7, MPI_DOUBLE_PRECISION, 0, 1 101, MPI_COMM_WORLD, status, ierr ) if (iam .ne. 0) call MPI_RECV(maxrow1,6, MPI_INTEGER, 0, 1 111, MPI_COMM_WORLD, status, ierr ) c dx2 = dx*dx dy2 = dy*dy c- c- Begining of time loop c- do 1000 itime = 1,maxtime c c- ctime = 0 itercnt = 0 done = 0 t = t + dt c if (iam .le. internod) then c- c- Region 1 c- ========----------------------------------------------- c- istef = 0 iallr = 0 c- c- if (iam .ne. 0) call MPI_RECV(temp1,transiz, 1 MPI_DOUBLE_PRECISION, 0, 1 130, MPI_COMM_WORLD, status, ierr ) do 100 i = 1,maxrow1/2 do 100 j = 1,maxcol1 u1t(i,j) = temp1(i,j) 100 continue c if (iam .ne. 0) call MPI_RECV(temp1,transiz, 1 MPI_DOUBLE_PRECISION, 0, 1 131, MPI_COMM_WORLD, status, ierr ) do 110 i = maxrow1/2 + 1, maxrow1 do 110 j = 1,maxcol1 u1t(i,j) = temp1(i,j) 110 continue c c- Initialize boundaries and initial guess for c- u1new,u1old c- c call MPI_BARRIER(MPI_COMM_WORLD, ierr) c- stime = MPI_WTIME() cstime = MPI_WTIME() c do 210 i = 2,maxrow1 do 220 j = 2,maxcol1 u1new(i,j) = u1t(i,j) u1old(i,j) = u1new(i,j) 220 continue 210 continue c- do 230 j = 1,maxcol1 u1new(1,j) = t u1old(1,j) = u1new(1,j) 230 continue c- do 240 i = 2,maxrow1 u1new(i,1) = t u1old(i,1) = u1new(i,1) 240 continue c- c- Begin the sor loop for solution at time t c- flag = 0 icount = 0 c- .......................................................... do while (done .eq. 0) itercnt = itercnt + 1 icount = icount +1 lastrow = (iam + 1)*numrows1 if (iam .eq. internod) then lastrow = lastrow - 1 + extra1 endif c- c- pass through top region up to and including bottom boundary c- if not the interface node. c- do 310 j = 2,maxcol1 - 1 do 300 i = iam*numrows1 + 2,lastrow u1new(i,j) = (1.0/r)*(1.0 + (u1t(i,j)/dt) + 1 (1.0/dx2)*(u1new(i,j+1) + u1new(i,j-1)) + 2 (1.0/dy2)*(u1new(i+1,j) + u1new(i-1,j))) u1new(i,j) = (1.0-w)*u1old(i,j) + w*u1new(i,j) 300 continue 310 continue c- c- equation for reflection of left hand side c- do 9300 i = iam*numrows1 + 2,lastrow u1new(i,maxcol1) = (1.0/r)*(1.0 + (u1t(i,maxcol1)/dt) + 1 (2.0/dx2)*u1new(i,maxcol1-1) + 2 (1.0/dy2)*(u1new(i+1,maxcol1) + u1new(i-1,maxcol1))) u1new(i,maxcol1) = (1.0-w)*u1old(i,maxcol1) + 1 w*u1new(i,maxcol1) c- 9300 continue c- c- Bottom of the region for interface node c- if (iam .eq. internod) then c- bottom of region up to a1-1 do 320 j = 2,a1-1 u1new(maxrow1,j) = (1.0/r)*(1.0 + (u1t(maxrow1,j)/dt) + 1 (1.0/dx2)*(u1new(maxrow1,j+1) + u1new(maxrow1,j-1)) + 2 (2.0/dy2)*u1new(maxrow1-1,j)) u1new(maxrow1,j) = (1.0-w)*u1old(maxrow1,j) + 1 w*u1new(maxrow1,j) 320 continue c- c- end of bottom region for internod c- endif c- do 350 j = 1 , maxcol1 tmp(j) = u1new(lastrow,j) 350 continue c- endif c- cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_SEND(tmp,vecsiz,MPI_DOUBLE_PRECISION,iam+1, 1 250, MPI_COMM_WORLD, ierr) cstime = MPI_WTIME() c- calculate interface for all nodes except first one c- if(iam .ne. 0) then cetime = MPI_WTIME() ctime = ctime + cetime -cstime call MPI_RECV(tmp,vecsiz,MPI_DOUBLE_PRECISION,iam-1, 1 250, MPI_COMM_WORLD, status, ierr) cstime = MPI_WTIME() i = iam*numrows1 + 1 do 360 j = 2,maxcol1 - 1 u1new(i,j) = (1.0/r)*(1.0 + (u1t(i,j)/dt) + 1 (1.0/dx2)*(u1new(i,j+1) + u1new(i,j-1)) + 2 (1.0/dy2)*(u1new(i+1,j) + tmp(j))) u1new(i,j) = (1.0-w)*u1old(i,j) + w*u1new(i,j) 360 continue c- c- equation for reflection of left hand side c- u1new(i,maxcol1) = (1.0/r)*(1.0 + (u1t(i,maxcol1)/dt) + 1 (2.0/dx2)*u1new(i,maxcol1-1) + 2 (1.0/dy2)*(u1new(i+1,maxcol1) + tmp(maxcol1))) u1new(i,maxcol1) = (1.0-w)*u1old(i,maxcol1) + 1 w*u1new(i,maxcol1) endif c- c- end of if loop for interface if not first node c- c- c- Send updated versions of u1new to previous iterate for next sor loop c- if ( iam .ne. 0) then do 400 j = 1,maxcol1 tmp(j) = u1new(iam*numrows1 + 1,j) 400 continue cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_SEND(tmp,vecsiz,MPI_DOUBLE_PRECISION,iam-1, 1 260, MPI_COMM_WORLD, ierr) cstime = MPI_WTIME() endif c- c- cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_RECV(tmp,vecsiz,MPI_DOUBLE_PRECISION,iam+1, 1 260, MPI_COMM_WORLD, status, ierr) cstime = MPI_WTIME() if (iam .ne. internod) then do 410 j = 1,maxcol1 u1new((iam+1)*numrows1+1,j) = tmp(j) 410 continue else do 420 l = a1,maxcol1 j = a2 + 4*(l-a1) u1new(maxrow1,l) = tmp(j) 420 continue lastrow = lastrow + 1 endif c- c- c- calculate error in the top region and update u1old err = 0.0 do 380 j = 2,maxcol1 do 370 i = iam*numrows1 + 1,lastrow if ( abs(u1new(i,j) - u1old(i,j)) .gt. err) then err = abs(u1new(i,j) - u1old(i,j)) endif u1old(i,j) = u1new(i,j) 370 continue 380 continue c- if (err.lt.tol) then flag = 1 endif c- cetime = MPI_WTIME() ctime = ctime + cetime - cstime c- if( mod(iallr,17).eq.0 ) then call MPI_ALLREDUCE(flag, temp, 1, MPI_INTEGER, 1 MPI_SUM, MPI_COMM_WORLD, ierr) endif iallr = iallr + 1 c- c cstime = MPI_WTIME() if ( temp .eq. maxnod + 1 ) then done = 1 endif c- c- if ( icount.gt.501) then done = 1 write(6,*) 'no convergency - nr. of iterations = ', 1 icount, ' message from node # ',iam stop endif c- c- flag = 0 c- enddo cetime = MPI_WTIME() etime = MPI_WTIME() time = etime - stime ctime = ctime + cetime - cstime if (iam .ne. 0) call MPI_SEND(time,3, 1 MPI_DOUBLE_PRECISION, 0, 1 400+iam, MPI_COMM_WORLD, ierr) if (iam .ne. 0) call MPI_SEND(itime,3,MPI_INTEGER, 0, 1 1400+iam, MPI_COMM_WORLD, ierr) c- do 500 j = 1,maxcol1 do 500 i = 1,maxrow1/2 temp1(i,j) = u1new(i,j) 500 continue if (iam .ne. 0) call MPI_SEND(temp1,transiz, 1 MPI_DOUBLE_PRECISION, 0, 1 300+iam, MPI_COMM_WORLD, ierr) c- do 510 j = 1,maxcol1 do 510 i = maxrow1/2 + 1,maxrow1 temp1(i,j) = u1new(i,j) 510 continue if (iam .ne. 0) call MPI_SEND(temp1,transiz, 1 MPI_DOUBLE_PRECISION, 0, 1 500+iam, MPI_COMM_WORLD, ierr) if (iam .eq. 0) return c- c- c- end of reg1 c- c-****************************************************************** c-****************************************************************** c- c- start of reg2 c- c- else c- c- if (iam .ne. 0) call MPI_RECV(u2t,transiz1, 1 MPI_DOUBLE_PRECISION, 0, 1 140, MPI_COMM_WORLD, status, ierr) c iallr = 0 d2x = dx/4.0 d2y = dy/4.0 d2x2 = d2x*d2x d2y2 = d2y*d2y r2 = 1.0/dt + 2.0/d2x2 + 2.0/d2y2 r3 = 1.0/dt + 2.0/d2x2 + 1.0/(2.0*d2y2) c- c- Initialize boundaries and initial guess for c- u2new,u2old c- c call MPI_BARRIER(MPI_COMM_WORLD, ierr) c stime = MPI_WTIME() cstime = MPI_WTIME() c do 1250 i = 1,maxrow2-1 do 1260 j = 2,maxcol2 u2new(i,j) = u2t(i,j) u2old(i,j) = u2new(i,j) 1260 continue 1250 continue c- do 1270 i = 1,maxrow2 u2new(i,1) = 0.0 u2old(i,1) = u2new(i,1) 1270 continue c- do 1280 j = 2,maxcol2 u2new(maxrow2,j) = 0.0 u2old(maxrow2,j) = u2new(maxrow2,j) 1280 continue c- c- Begin the sor loop for solution at time t c- flag = 0 icount = 0 do while (done .eq. 0) c- .................................................... c- itercnt = itercnt + 1 icount = icount + 1 lastrow = (iam - internod)*numrows2 c if ( iam .eq. maxnod) then lastrow = lastrow - 1 + extra2 endif c- c- c- calculate the bottom region c- c- do 1300 i = (iam - internod - 1)*numrows2 + 2,lastrow do 1310 j = 2,maxcol2-1 do 1300 i = (iam - internod - 1)*numrows2 + 2,lastrow u2new(i,j) = (1.0/r2)*(-b + (u2t(i,j)/dt) + 1 (1.0/d2x2)*(u2new(i,j+1) + u2new(i,j-1)) + 2 (1.0/d2y2)*(u2new(i+1,j) + u2new(i-1,j))) u2new(i,j) = (1.0-w)*u2old(i,j) + w*u2new(i,j) if (u2new(i,j) .lt. 0.0) then u2new(i,j) = 0.0 endif 1300 continue 1310 continue c - do 7300 i = (iam - internod - 1)*numrows2 + 2,lastrow u2new(i,maxcol2) = (1.0/r2)*(-b + (u2t(i,maxcol2)/dt) + 1 (2.0/d2x2)*u2new(i,maxcol2-1) + 2 (1.0/d2y2)*(u2new(i+1,maxcol2) + u2new(i-1,maxcol2))) u2new(i,maxcol2) = (1.0-w)*u2old(i,maxcol2) + 1 w*u2new(i,maxcol2) if (u2new(i,maxcol2) .lt. 0.0) then u2new(i,maxcol2) = 0.0 endif 7300 continue c- c- c- passing data along the interface c- if(iam .ne. maxnod) then do 1320 j = 1,maxcol2 tmp(j) = u2new(lastrow,j) 1320 continue cetime = MPI_WTIME() ctime = ctime + cetime - cstime c- c- c- call MPI_SEND(tmp, vecsiz, MPI_DOUBLE_PRECISION, iam+1, 1 250, MPI_COMM_WORLD, ierr) cstime = MPI_WTIME() endif c- c- receive information along interfaces c- cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_RECV(tmp, vecsiz, MPI_DOUBLE_PRECISION, iam-1, 1 250, MPI_COMM_WORLD, status, ierr) cstime = MPI_WTIME() c- c- .................................... if (iam .ne. internod+1) then c- c- i = (iam - internod -1)*numrows2 + 1 do 1330 j = 2,maxcol2-1 u2new(i,j) = (1.0/r2)*(-b + (u2t(i,j)/dt) + 1 (1.0/d2x2)*(u2new(i,j+1) + u2new(i,j-1)) + 2 (1.0/d2y2)*(u2new(i+1,j) + tmp(j))) u2new(i,j) = (1.0-w)*u2old(i,j) + w*u2new(i,j) if (u2new(i,j) .lt. 0.0) then u2new(i,j) = 0.0 endif 1330 continue u2new(i,maxcol2) = (1.0/r2)*(-b + (u2t(i,maxcol2)/dt) + 1 (2.0/d2x2)*u2new(i,maxcol2-1) + 2 (1.0/d2y2)*(u2new(i+1,maxcol2) + tmp(maxcol2))) u2new(i,maxcol2) = (1.0-w)*u2old(i,maxcol2) + 1 w*u2new(i,maxcol2) if (u2new(i,maxcol2) .lt. 0.0) then u2new(i,maxcol2) = 0.0 endif c- c- c- .................... else c- .................... c- interface of region from maxcol2 to 1 c- c- c- calculate information and interpolate along wet region c- u2new(1,a2) = 0.5*(u2new(5,a2) + tmp(a1)) if ( u2new(1,a2) .lt. 0.0 ) then u2new(1,a2) = 0.0 endif c- c- c- do 340 l = a1+1, maxcol1 j = a2 + 4*(l - a1) u2new(1,j) = 0.5*(u2new(5,j) + tmp(l)) if (u2new(1,j) .lt. 0.0 ) then u2new(1,j) = 0.0 endif c- c- interpolate middle nodes c- ind2 = j - 4 deltay = 0.25*(u2new(1,j) - u2new(1,j-4)) do 325 i = 1,3 ind2 = ind2 + 1 u2new(1,ind2) = u2new(1,j-4) + deltay if (u2new(1,ind2) .lt. 0.0) then u2new(1,ind2) = 0.0 endif 325 continue 340 continue c- c- c- calculate equation for u along bottom of mask c- do 1350 j = a2-1,2,-1 u2new(1,j) = (1.0/r2)*(-b + (u2t(1,j)/dt) + 1 (1.0/d2x2)*(u2new(1,j+1) + u2new(1,j-1)) + 2 (2.0/d2y2)*u2new(2,j)) u2new(1,j) = (1.0-w)*u2old(1,j) + w*u2new(1,j) if (u2new(1,j) .lt. 0.0) then u2new(1,j) = 0.0 endif 1350 continue c- endif c- ............................................... c- c- Update u2new, and u1new from the node below for the last + 1 row c- c- row = (iam-internod-1)*numrows2 + 1 do 1420 j = 1,maxcol2 tmp(j) = u2new(row,j) 1420 continue cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_SEND(tmp, vecsiz, MPI_DOUBLE_PRECISION, iam-1, 1 260, MPI_COMM_WORLD, ierr) cstime = MPI_WTIME() c- c- if(iam.ne.maxnod) then cetime = MPI_WTIME() ctime = ctime + cetime - cstime call MPI_RECV(tmp, vecsiz, MPI_DOUBLE_PRECISION, iam+1, 1 260, MPI_COMM_WORLD, status, ierr) cstime = MPI_WTIME() do 1430 j = 1,maxcol2 u2new(lastrow+1,j) = tmp(j) 1430 continue endif c- c- c- calculate error in bottom region and update u2old c- err = 0.0 do 1410 j = 2,maxcol2 do 1400 i = (iam-internod-1)*numrows2 +1,lastrow if ( abs(u2new(i,j) - u2old(i,j)) .gt. err) then err = abs(u2new(i,j) - u2old(i,j)) endif u2old(i,j) = u2new(i,j) 1400 continue 1410 continue c- if(err .lt. tol) then flag = 1 endif c- c- c- sychronize all nodes to see if all have converged c- c- cetime = MPI_WTIME() ctime = ctime + cetime - cstime c- if( mod(iallr,17).eq.0) then call MPI_ALLREDUCE(flag, temp, 1, MPI_INTEGER, 1 MPI_SUM, MPI_COMM_WORLD, ierr) endif iallr = iallr + 1 c- c cstime = MPI_WTIME() if(temp .eq. maxnod+1) then done = 1 endif c- c- if ( icount.gt.501) then done = 1 write(6,*) 'no convergency - nr. of iterations = ', 1 icount, ' message from node # ',iam stop endif c- c- flag = 0 c- c- end of sor loop c- enddo cetime = MPI_WTIME() etime = MPI_WTIME() ctime = ctime + cetime - cstime time = etime - stime call MPI_SEND(itime, 3, MPI_INTEGER, 0, 1 1400+iam, MPI_COMM_WORLD, ierr) call MPI_SEND(time, 3, 1 MPI_DOUBLE_PRECISION, 0, 1 400+iam, MPI_COMM_WORLD, ierr) call MPI_SEND(u2new, transiz1, MPI_DOUBLE_PRECISION, 0, 1 300+iam, MPI_COMM_WORLD, ierr) c- c- c- end of if loop for which region c- c- endif c c if(itime.gt.63) then c write(6,*) 'message from node#',iam c 1 , '--- end loop ---------time step reg.1 = ', itime c endif c- 1000 continue c write(6,*) ' now at the end of etch6n ' return end