diff --git a/m_compact.f90 b/m_compact.f90 index 4d7cfe6..9302cd6 100644 --- a/m_compact.f90 +++ b/m_compact.f90 @@ -8,10 +8,13 @@ lyf,lys,wyf,wys, & lzf,lzs,wzf,wzs ! lyzf,lyzs,wyzf,wyzs + + REAL*8, DIMENSION(:), ALLOCATABLE :: zz1, zz2, vz1, vz2 + INTEGER :: nxc,nyc,nzc REAL*8, PARAMETER :: ezero = 1.0e-14 - PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp + PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp,pdfp,pd2fp CONTAINS @@ -68,6 +71,22 @@ call nonp_lud(3,nzc) ENDIF + ALLOCATE(vz1(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + + ALLOCATE(zz1(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + + ALLOCATE(vz2(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + + ALLOCATE(zz2(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + + call ctdlu (3.0, nzc, vz1, zz1) + + call ctdlu (11.0/2.0, nzc, vz2, zz2) + ! CALL x_lud ! CALL yz_lud @@ -144,6 +163,232 @@ l(n)=1./d END SUBROUTINE ptdlu + SUBROUTINE ctdlu(lambda,n,v,z) + INTEGER, INTENT(IN) :: n + REAL*8, INTENT(IN) :: lambda + REAL*8, INTENT(OUT) :: v(n) + REAL*8, INTENT(OUT) :: z(n) + INTEGER :: i + REAL*8 :: sigma, alpha + + v(1)=0.0 + v(2)=lambda + DO i=3,n + v(i)=lambda - 1./v(i-1) + ENDDO + + if (lambda > 0.0) then + alpha = (-lambda + sqrt(lambda**2 - 4.0))/2.0 + else + alpha = (-lambda - sqrt(lambda**2 - 4.0))/2.0 + endif + + sigma = (1.0 + alpha**2) / (lambda * (1.0 - alpha**2) * (1.0 - alpha**n)) + + DO i=1,n + z(i) = sigma * (alpha**(i-1) + alpha**(n-i+1)) + ENDDO + END SUBROUTINE ctdlu + + + SUBROUTINE pdfp(x,xu,xl,dx,h,nd,n,nall,dir) + + INTEGER,INTENT(IN) :: nall,n,nd,dir + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,2) :: xu, xl + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx + + INTEGER :: i,j + REAL*8 :: r1,r2,h1 + + h1=1./h + r1=7./3. + r2=1./12. + + !if (myid.eq.0) write(*,*) "parallel dfp commuication" + if (myid.eq.master) then + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + elseif (myid.eq.numprocs-1) then + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + else + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + endif + + DO j=1,nd + dx(j,1) =(r1*( x(j,2)-xl(j,2)) +r2*( x(j,3)-xl(j,1))) *h1 + ENDDO + + DO j=1,nd + dx(j,2) =(r1*( x(j,3)- x(j,1)) +r2*( x(j,4)-xl(j,2))) *h1 + ENDDO + + DO j=1,nd + dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1 + ENDDO + + DO j=1,nd + dx(j,n) =(r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1 + ENDDO + + DO i=3,n-2 + DO j=1,nd + dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))*h1 + ENDDO + ENDDO + + + IF (dir.eq.1) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction + IF (dir.eq.2) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction + IF (dir.eq.3) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction + + END SUBROUTINE pdfp + + + SUBROUTINE pd2fp(x,xu,xl,dx,h,nd,n,nall,dir) + INTEGER,INTENT(IN) :: nall,n,nd,dir + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL*8 :: h2,r1,r2,t1,t2 + + REAL*8,DIMENSION(nd,2) :: xu, xl + + h2=1./(h*h) + r1=6. + r2=3./8. + + !if (myid.eq.0) write(*,*) "parallel d2fp commuication" + if (myid.eq.master) then + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, numprocs-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, numprocs-1, numprocs-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + elseif (myid.eq.numprocs-1) then + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, 0, 0, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + else + + call MPI_ISEND (x(1,1), 2*nd, MPI_REAL8, myid-1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xu, 2*nd, MPI_REAL8, myid+1, myid+1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, & + MPI_COMM_TASK, mpi_request, mpi_err) + + call MPI_RECV (xl, 2*nd, MPI_REAL8, myid-1, myid-1, & + MPI_COMM_TASK, mpi_status, mpi_err) + + endif + + + DO j=1,nd + + t1 = (x(j,2)-2.*x(j,1)+xl(j,2)) + t2 = (x(j,3)-2.*x(j,1)+xl(j,1)) + dx(j,1)=(r1*t1+r2*t2)*h2 + + ENDDO + + DO j=1,nd + + t1 = (x(j,3)-2.*x(j,2)+ x(j,1)) + t2 = (x(j,4)-2.*x(j,2)+xl(j,2)) + dx(j,2)=(r1*t1+r2*t2)*h2 + + ENDDO + + DO j=1,nd + + t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2)) + t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3)) + dx(j,n-1)=(r1*t1+r2*t2)*h2 + + ENDDO + + DO j=1,nd + + t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1)) + t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2)) + dx(j,n)=(r1*t1+r2*t2)*h2 + + ENDDO + + DO i=3,n-2 + DO j=1,nd + + t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1)) + t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2)) + dx(j,i)=(r1*t1+r2*t2)*h2 + + ENDDO + ENDDO + + IF (dir.eq.1) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! x-direction + IF (dir.eq.2) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! y-direction + IF (dir.eq.3) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! z-direction + END SUBROUTINE pd2fp + + SUBROUTINE dfnonp(n,h,x,dx,nd,dir) INTEGER,INTENT(IN) :: n,nd,dir REAL*8,INTENT(IN) :: h @@ -329,6 +574,187 @@ ENDDO END SUBROUTINE ptdslv + SUBROUTINE ctdslv(r,va,za,nd,n,nall) + INTEGER,PARAMETER :: nb = 256 + INTEGER,INTENT(IN) :: n,nd,nall + REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL*8,INTENT(IN),DIMENSION(nall) :: va + REAL*8,INTENT(IN),DIMENSION(nall) :: za + REAL*8,DIMENSION(n) :: v,z + INTEGER i,j + INTEGER ii,jj + REAL*8, DIMENSION(nd) :: sum + + REAL*8, DIMENSION(nd) :: x1, x1_ + + INTEGER npart, nbase, nlow, nupp + INTEGER :: pid, np + REAL*8, DIMENSION(nb) :: r0, r1, sum0 + REAL*8, DIMENSION(nb,2) :: buf + + pid = myid + np = numprocs + + npart = nall / np + nbase = pid * npart + nlow = pid * npart + 1 + nupp = (pid + 1) * npart + + v = va(nlow:nupp) + z = za(nlow:nupp) + + if (npart.lt.4) then + ! assertion fail + endif + + x1 = 0.0 + x1_ = 0.0 + + do i=1,n + x1_ = x1_ + z(i)*r(:,i) + enddo + + ! call MPI_ALLREDUCE(x1_, x1, nd, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err) + call MPI_REDUCE(x1_, x1, nd, MPI_REAL8, MPI_SUM, master, MPI_COMM_TASK, mpi_err) + + ! first process + if (myid.eq.0) then + + CALL MPI_ISEND(x1, nd, MPI_REAL8, np-1, 100, MPI_COMM_TASK, mpi_request, mpi_err) + + r(:,1) = x1 + r(:,2) = r(:,2) - x1 + + ! last process + elseif (pid.eq.(np-1)) then + + CALL MPI_RECV(x1, nd, MPI_REAL8, master, 100, MPI_COMM_TASK, mpi_status, mpi_err) + + r(:,n) = r(:,n) - x1 + + ! intermediate process + else + + endif + + ! first process + if (myid.eq.0) then + + DO jj=1,nd,nb + + DO j=jj,jj+nb-1 + r(j,2)=r(j,2)/v(2) + ENDDO + + DO i=3,n + DO j=jj,jj+nb-1 + r(j,i)=(r(j,i)-r(j,i-1))/v(i) + ENDDO + ENDDO + + r0 = r(jj:jj+nb-1,n) + CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err) + + ENDDO + + DO jj=1,nd,nb + + CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) + + i=n + DO j=jj,jj+nb-1 + r(j,i) = r(j,i) - r1(j-jj+1)/v(i) + ENDDO + + DO i=n-1,2,-1 + DO j=jj,jj+nb-1 + r(j,i) = r(j,i) - r(j,i+1)/v(i) + ENDDO + ENDDO + + ENDDO + + ! last process + elseif (pid.eq.(np-1)) then + + DO jj=1,nd,nb + + CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + + i=1 + DO j=jj,jj+nb-1 + r(j,i)=(r(j,i) - r1(j-jj+1))/v(i) + ENDDO + + DO i=2,n + DO j=jj,jj+nb-1 + r(j,i)=(r(j,i) - r(j,i-1))/v(i) + ENDDO + ENDDO + + ENDDO + + DO jj=1,nd,nb + + DO i=n-1,1,-1 + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-r(j,i+1)/v(i) + ENDDO + ENDDO + + r0 = r(jj:jj+nb-1,1) + CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err) + + ENDDO + + ! intermediate process + else + + DO jj=1,nd,nb + + CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err) + + i=1 + DO j=jj,jj+nb-1 + r(j,i)=(r(j,i)-r1(j-jj+1))/v(i) + ENDDO + + DO i=2,n + DO j=jj,jj+nb-1 + r(j,i)=(r(j,i)-r(j,i-1))/v(i) + ENDDO + ENDDO + + r0 = r(jj:jj+nb-1,n) + CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err) + + ENDDO + + DO jj=1,nd,nb + + CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err) + + i=n + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-r1(j-jj+1)/v(i) + ENDDO + + DO i=n-1,1,-1 + DO j=jj,jj+nb-1 + r(j,i)=r(j,i)-r(j,i+1)/v(i) + ENDDO + ENDDO + + r0 = r(jj:jj+nb-1,1) + CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err) + + ENDDO + + endif + + END SUBROUTINE ctdslv + + SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall) INTEGER,PARAMETER :: nb = 64 INTEGER,INTENT(IN) :: n,nd,nall diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 375ca44..0ce9382 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -546,8 +546,8 @@ module m_fdm_calc ENDDO ENDDO - CALL par_dfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3) - CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3) + CALL pdfp (fz, fzl, fzu, dfz, hy, 4*xx*yy, zz, yy, 3) + CALL pd2fp(fzz, fzzl, fzzu, dfzz, hy, xx*yy, zz, yy, 3) DO k=1,zz DO j=1,yy