From d947ffe62a86740a58937012a66d1d2b3513b050 Mon Sep 17 00:00:00 2001 From: ignis Date: Fri, 5 Mar 2021 22:45:26 +0900 Subject: [PATCH] added extract ddxc, timestamp to log numpy memmap py3 compatibility fix --- Compact.f90 | 539 ++++++++++++++++++++++++++++++++++++++++++++++++ compact_doc.py | 2 + extract_c.py | 15 +- extract_ddxc.py | 33 +++ extract_u.py | 15 +- pycompact.py | 528 +++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 1120 insertions(+), 12 deletions(-) create mode 100644 Compact.f90 create mode 100644 compact_doc.py create mode 100644 extract_ddxc.py create mode 100644 pycompact.py diff --git a/Compact.f90 b/Compact.f90 new file mode 100644 index 0000000..d2ee2b7 --- /dev/null +++ b/Compact.f90 @@ -0,0 +1,539 @@ + MODULE Compact + IMPLICIT NONE + + REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & + lyf,lys,wyf,wys, & + lzf,lzs,wzf,wzs + + INTEGER :: nxc,nyc,nzc + + REAL(KIND=8), PARAMETER :: ezero = 1.0e-14 + + CONTAINS + + SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz + + CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp) + + CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp) + + END SUBROUTINE ludcmp + + SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz + +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + ALLOCATE(lxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(xp.eq.0) THEN + ALLOCATE(wxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ENDIF + + ALLOCATE(lyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(yp.eq.0) THEN + ALLOCATE(wyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ENDIF + + ALLOCATE(lzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(zp.eq.0) THEN + ALLOCATE(wzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ENDIF + + END SUBROUTINE ludcmp_allocate + + SUBROUTINE ludcmp_deallocate(xp,yp,zp) + INTEGER, INTENT(IN) :: xp,yp,zp + +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + DEALLOCATE(lxf) + DEALLOCATE(lxs) + IF(xp.eq.0) THEN + DEALLOCATE(wxf) + DEALLOCATE(wxs) + ENDIF + + DEALLOCATE(lyf) + DEALLOCATE(lys) + IF(yp.eq.0) THEN + DEALLOCATE(wyf) + DEALLOCATE(wys) + ENDIF + + DEALLOCATE(lzf) + DEALLOCATE(lzs) + IF(zp.eq.0) THEN + DEALLOCATE(wzf) + DEALLOCATE(wzs) + ENDIF + + END SUBROUTINE ludcmp_deallocate + + SUBROUTINE ludcmp_testalloc + + IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated" + IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated" + IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated" + IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated" + + IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated" + IF (.not. ALLOCATED(lys)) print *, "lys not allocated" + IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated" + IF (.not. ALLOCATED(wys)) print *, "wys not allocated" + + IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated" + IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated" + IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated" + IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated" + + END SUBROUTINE ludcmp_testalloc + + SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz + + ! CALL ludcmp_testalloc + +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + IF(xp.eq.0) THEN + CALL p_lud(1,nxc) + ELSE + CALL nonp_lud(1,nxc) + ENDIF + + IF(yp.eq.0) THEN + CALL p_lud(2,nyc) + ELSE + call nonp_lud(2,nyc) + ENDIF + + IF(zp.eq.0) THEN + CALL p_lud(3,nzc) + ELSE + call nonp_lud(3,nzc) + ENDIF + + END SUBROUTINE ludcmp_calculate + + + SUBROUTINE test_nonp_lud1(xx, coef) + INTEGER :: xx + REAL(KIND=8), DIMENSION(xx) :: aa + REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef + aa=3. + aa(1)=0.5 ; aa(2)=4. + aa(xx-1)=4. ; aa(xx)=0.5 + + CALL stdlu(aa,xx,coef) + END SUBROUTINE test_nonp_lud1 + + + SUBROUTINE test_nonp_lud2(xx, coef) + INTEGER :: xx + REAL(KIND=8), DIMENSION(xx) :: aa + REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef + aa=5.5 + aa(1)=2./11. ; aa(2)=10. + aa(xx-1)=10. ; aa(xx)=2./11. + + CALL stdlu(aa,xx,coef) + END SUBROUTINE test_nonp_lud2 + + + SUBROUTINE test_p_lud1(xx, coef1, coef2) + INTEGER :: xx + REAL(KIND=8) :: a + REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2 + a=3. ! first derivative + + CALL ptdlu(a,xx,coef1,coef2) ! x-direction + END SUBROUTINE test_p_lud1 + + + SUBROUTINE test_p_lud2(xx, coef1, coef2) + INTEGER :: xx + REAL(KIND=8) :: a + REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2 + a=11./2. ! second derivative + + CALL ptdlu(a,xx,coef1,coef2) ! x-direction + END SUBROUTINE test_p_lud2 + + + SUBROUTINE nonp_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL(KIND=8), DIMENSION(xx) :: aa + aa=3. + aa(1)=0.5 ; aa(2)=4. + aa(xx-1)=4. ; aa(xx)=0.5 +! first derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction + aa=5.5 + aa(1)=2./11. ; aa(2)=10. + aa(xx-1)=10. ; aa(xx)=2./11. +! second derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction + END SUBROUTINE nonp_lud + + + SUBROUTINE p_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL(KIND=8) :: a + a=3. ! first derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction + a=11./2. ! second derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction + END SUBROUTINE p_lud + + SUBROUTINE stdlu(a,n,l) + INTEGER :: n + REAL(KIND=8), INTENT(IN) :: a(n) + REAL(KIND=8), INTENT(OUT) :: l(n) + REAL(KIND=8) :: d + INTEGER :: i + l(1)=1.0d0/a(1) + DO i=2,n + d=a(i)-l(i-1) + l(i)=1.0d0/d + ENDDO + END SUBROUTINE stdlu + + SUBROUTINE ptdlu(a,n,l,w) + INTEGER :: n + REAL(KIND=8), INTENT(OUT) :: a + REAL(KIND=8), INTENT(OUT) :: l(n),w(n) + INTEGER :: i + REAL(KIND=8) :: aa(n),d + + DO i=1,n-1 + aa(i)=a + ENDDO + i=n-1 + call stdlu(aa,i,l) + w(1)=1.0 + DO i=2,n-2 + w(i)=-l(i-1)*w(i-1) + ENDDO + w(n-1)=1.0-l(n-2)*w(n-2) + DO i=1,n-1 + w(i)=w(i)*l(i) + ENDDO + d=a + DO i=1,n-1 + d=d-w(i)*w(i)/l(i) + ENDDO + l(n)=1./d + END SUBROUTINE ptdlu + + + SUBROUTINE rhs1np(n,h,x,dx,nd) + INTEGER,INTENT(IN) :: n,nd + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + + h1=1.d0/h + + r1=7.d0/3.d0 + r2=1.d0/12.d0 + r3=3. + a=-1.25 + b=1. + c=0.25 + + DO j=1,nd + dx(j,n-1)=x(j,n)-x(j,n-2) + dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)) + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)) + dx(j,2)=x(j,3)-x(j,1) + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0. + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0. + dx(j,n-1)=dx(j,n-1)*h1*r3 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(j,2)*h1*r3 + ENDDO + + DO i=3,n-2 + DO j=1,nd + t1=x(j,i+1)-x(j,i-1) + t2=x(j,i+2)-x(j,i-2) + dx(j,i)=h1*(r1*t1+r2*t2) + ENDDO + ENDDO + + END SUBROUTINE rhs1np + + + SUBROUTINE dfnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + + CALL rhs1np (n,h,x,dx,nd) + + IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction + END SUBROUTINE dfnonp + + SUBROUTINE dfp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: r1,r2,h1 + + + ! print *, "dfnonp received (nd,n)", nd, n + + h1=1./h + r1=7./3. + r2=1./12. + + DO j=1,nd + dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3))) + dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2))) + dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1))) + dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n))) + dx(j,n-1)=dx(j,n-1)*h1 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(j,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))) + dx(j,i)=dx(j,i)*h1 + ENDDO + ENDDO + + IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction + + END SUBROUTINE dfp + + SUBROUTINE ptdslv(r,n,l,w,nd) + INTEGER,INTENT(IN) :: n,nd + REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r + REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w + INTEGER i,j + REAL(KIND=8), DIMENSION(nd) :: sum + DO j=1,nd + sum(j)=w(1)*r(j,1) + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n-1 + DO j=1,nd + r(j,i)=r(j,i)-r(j,i-1) + sum(j)=sum(j)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO + ENDDO + DO j=1,nd + r(j,n)=l(n)*(r(j,n)-sum(j)) + r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n) + ENDDO + DO i=n-2,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n) + ENDDO + ENDDO + END SUBROUTINE ptdslv + + SUBROUTINE d2fp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: h2,r1,r2,t1,t2 + + + h2=1./(h*h) + r1=6. + r2=3./8. + DO j=1,nd + t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2)) + t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3)) + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0. + IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0. + dx(j,n-1)=(r1*t1+r2*t2) + + t1 = (x(j,1)-2.*x(j,n)+x(j,n-1)) + t2 = (x(j,2)-2.*x(j,n)+x(j,n-2)) + IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0. + IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0. +! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) & +! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2))) + dx(j,n)=(r1*t1+r2*t2) + + t1 = (x(j,2)-2.*x(j,1)+x(j,n)) + t2 = (x(j,3)-2.*x(j,1)+x(j,n-1)) + IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0. + IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0. +! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) & +! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1))) + dx(j,1)=(r1*t1+r2*t2) + + t1 = (x(j,3)-2.*x(j,2)+x(j,1)) + t2 = (x(j,4)-2.*x(j,2)+x(j,n)) + IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0. + IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0. +! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) & +! +r2*(x(j,4)-2.*x(j,2)+x(j,n))) + dx(j,2)=(r1*t1+r2*t2) + + dx(j,n-1)=dx(j,n-1)*h2 + dx(j,n)=dx(j,n)*h2 + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*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)) + IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0. + IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0. +! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) & +! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2))) + dx(j,i)=(r1*t1+r2*t2) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction + END SUBROUTINE d2fp + + SUBROUTINE tdslv(r,n,l,nd) + INTEGER,INTENT(IN) :: n,nd + REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r + REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l + INTEGER i,j + REAL(KIND=8) t1 + DO j=1,nd + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n + DO j=1,nd + t1=r(j,i)-r(j,i-1) + r(j,i)=l(i)*t1 + ENDDO + ENDDO + DO i=n-1,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1) + ENDDO + ENDDO + END SUBROUTINE tdslv + + SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL(KIND=8),INTENT(IN) :: h + REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x + REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2 + + + h2=1./(h*h) + r1=6. + r2=3./8. + r3=12. + a=13./11. + b=-27./11. + c=15./11. + e=-1./11. + + DO j=1,nd + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4)) + dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1)) + dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2)) + dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3)) + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0. + IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0. + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*h2*r3 + dx(j,n-1)=dx(j,n-1)*h2*r3 + dx(j,n)=dx(j,n)*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)) + IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0. + IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0. + +! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) & +! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2))) + dx(j,i)=(r1*t1+r2*t2) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + + IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction + + END SUBROUTINE d2fnonp + + END MODULE Compact diff --git a/compact_doc.py b/compact_doc.py new file mode 100644 index 0000000..f25cb7c --- /dev/null +++ b/compact_doc.py @@ -0,0 +1,2 @@ +from compact import compact +print compact.__doc__ diff --git a/extract_c.py b/extract_c.py index 7dfac11..9fcd1c5 100644 --- a/extract_c.py +++ b/extract_c.py @@ -1,4 +1,5 @@ import argparse +import datetime import numpy as np import dnstool @@ -17,11 +18,13 @@ case = cases[casename] nx, ny, nz = case.shape -with open("c.dat", "w+") as ufile: +ufile = "c.dat" - storage = np.memmap(ufile, dtype=np.double, shape=(nx, len(case.data_files), ny, nz)) +storage = np.memmap(ufile, mode='w+', dtype=np.double, shape=(nx, len(case.data_files), ny, nz)) - for i, fname in enumerate(case.data_files): - print(i, fname) - time, y = case.read_single_field(fname, field_idx) - storage[:,i,:,:] = 1. - y.T +for i, fname in enumerate(case.data_files): + print(datetime.datetime.now(), i, fname) + time, y = case.read_single_field(fname, field_idx) + storage[:,i,:,:] = 1. - y.T + +storage.flush() diff --git a/extract_ddxc.py b/extract_ddxc.py new file mode 100644 index 0000000..86bab84 --- /dev/null +++ b/extract_ddxc.py @@ -0,0 +1,33 @@ +import argparse +import datetime +import numpy as np +import dnstool +from pycompact import CompactScheme + +# Commandline argument parser +parser = argparse.ArgumentParser() +parser.add_argument("-c", "--case", help="target case name", required=True) +args = parser.parse_args() +params = vars(args) + +casename = params["case"] +field_idx = 4 + +cases = dnstool.case_library() + +case = cases[casename] + +nx, ny, nz = case.shape + +cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2) + +ufile = "ddxc.dat" + +storage = np.memmap(ufile, mode='w+', dtype=np.double, shape=(nx, len(case.data_files), ny, nz)) + +for i, fname in enumerate(case.data_files): + print(datetime.datetime.now(), i, fname) + time, y = case.read_single_field(fname, field_idx) + storage[:,i,:,:] = cs.ddx(1. - y).T + +storage.flush() diff --git a/extract_u.py b/extract_u.py index d25f0a8..69d72b4 100644 --- a/extract_u.py +++ b/extract_u.py @@ -1,4 +1,5 @@ import argparse +import datetime import numpy as np import dnstool @@ -17,11 +18,13 @@ case = cases[casename] nx, ny, nz = case.shape -with open("u.dat", "w+") as ufile: +ufile = "u.dat" - storage = np.memmap(ufile, dtype=np.double, shape=(nx, len(case.data_files), ny, nz)) +storage = np.memmap(ufile, mode='w+', dtype=np.double, shape=(nx, len(case.data_files), ny, nz)) - for i, fname in enumerate(case.data_files): - print(i, fname) - time, u = case.read_single_field(fname, field_idx) - storage[:,i,:,:] = u.T +for i, fname in enumerate(case.data_files): + print(datetime.datetime.now(), i, fname) + time, u = case.read_single_field(fname, field_idx) + storage[:,i,:,:] = u.T + +storage.flush() diff --git a/pycompact.py b/pycompact.py new file mode 100644 index 0000000..0b46dab --- /dev/null +++ b/pycompact.py @@ -0,0 +1,528 @@ +import numpy as np +from compact import compact + +class CompactScheme: + + def __init__ (self, nx, ny, nz, px, py, pz, lx, ly, lz): + + pi8 = np.arccos(-1., dtype=np.float64) + + self.shape = (nz, ny, nx) + + self.px = px + self.py = py + self.pz = pz + + h = pi8 * lx / nx + self.hx = h + self.hy = h + self.hz = h + + # Allocate LU + compact.lxf = np.zeros(nx, dtype=np.float64) + compact.lxs = np.zeros(nx, dtype=np.float64) + compact.wxf = np.zeros(nx, dtype=np.float64) + compact.wxs = np.zeros(nx, dtype=np.float64) + + compact.lyf = np.zeros(ny, dtype=np.float64) + compact.lys = np.zeros(ny, dtype=np.float64) + compact.wyf = np.zeros(ny, dtype=np.float64) + compact.wys = np.zeros(ny, dtype=np.float64) + + compact.lzf = np.zeros(nz, dtype=np.float64) + compact.lzs = np.zeros(nz, dtype=np.float64) + compact.wzf = np.zeros(nz, dtype=np.float64) + compact.wzs = np.zeros(nz, dtype=np.float64) + + bcx = 0 if px else 1 + bcy = 0 if py else 1 + bcz = 0 if pz else 1 + + compact.ludcmp_calculate(nx, ny, nz, bcx, bcy, bcz) + + + def test_ludcmp (self): + import pprint + pp = pprint.PrettyPrinter(indent=4) + + # First Derivative Non-periodic BC + l1 = compact.test_nonp_lud1(self.shape[-1]) + + # Second Derivative Non-periodic BC + l2 = compact.test_nonp_lud2(self.shape[-1]) + + print ("Test Internally Calculated Non-periodic Coefs") + print (np.linalg.norm((l1 - compact.lxf)/compact.lxf)) + print (np.linalg.norm((l2 - compact.lxs)/compact.lxs)) + + + def py_rhs_1_np (self, x): + dx = np.zeros(x.shape) + + h1 = 1./self.hx + + r1 = 7./3. + r2 = 1./12. + r3 = 3. + a = -1.25 + b = 1. + c = 0.25 + + nd, n = x.shape + + dx[:, -2] = x[:, -1] - x[:, -3] + dx[:, -1] = - (a*x[:, -1] + b*x[:, -2] + c*x[:, -3]) + dx[:, 0] = (a*x[:, 0] + b*x[:, 1] + c*x[:, 2]) + dx[:, 1] = x[:, 2] - x[:, 0] + + dx[:,-2] = dx[:,-2]*h1*r3 + dx[:,-1] = dx[:,-1]*h1 + dx[:,0] = dx[:,0]*h1 + dx[:,1] = dx[:,1]*h1*r3 + + for i in range(2,n-2): + t1=x[:,i+1]-x[:,i-1] + t2=x[:,i+2]-x[:,i-2] + dx[:,i]=h1*(r1*t1+r2*t2) + + return dx + + + def py_tdslv(self, r, l): + nd, n = r.shape + + r[:,0] = r[:,0] * l[0] + + for i in range(1,n): + r[:,i] = l[i] * (r[:,i] - r[:,i-1]) + + for i in range(n-1)[::-1]: + r[:,i] = r[:,i] - l[i] * r[:,i+1] + + + def test_dfnonp (self): + + x = np.sin(1.1 * np.arange(512) * self.hx).reshape((1,-1)) + + exact = 1.1 * np.cos(1.1 * np.arange(512) * self.hx).reshape((1,-1)) + + + + print ("First Non-periodic RHS Test") + + dx = self.py_rhs_1_np(x) + + dx_fortran = compact.rhs1np(self.hx, x) + + print ("RelError Norm: ", np.linalg.norm((dx - dx_fortran) / dx_fortran)) + print ("RelError Min : ", ((dx - dx_fortran) / dx_fortran).min()) + print ("RelError Min : ", ((dx - dx_fortran) / dx_fortran).max()) + + + + + print ("First Non-periodic TD SOLVE Test") + + l1 = compact.test_nonp_lud1(512) + + self.py_tdslv(dx, l1) + + compact.tdslv(dx_fortran,l1) + + print ("dx - exact") + print ("RelError Norm: ", np.linalg.norm((dx - exact) / exact)) + print ("RelError Min : ", ((dx - exact) / exact).min()) + print ("RelError Min : ", ((dx - exact) / exact).max()) + + print ("dx_fortran - exact") + print ("RelError Norm: ", np.linalg.norm((dx_fortran - exact) / exact)) + print ("RelError Min : ", ((dx_fortran - exact) / exact).min()) + print ("RelError Min : ", ((dx_fortran - exact) / exact).max()) + + ''' + import pprint + pp = pprint.PrettyPrinter(indent=4) + + pp.pprint ((dx - exact) / exact) + + pp.pprint ((zip ((dx - exact).ravel(), dx.ravel(), exact.ravel()))) + ''' + + def verify_nonp_lud1(self): + + print ("Non-periodic coef first derivative") + + nx = 512 + aa = np.ones(nx) * 3. + aa[0] = 0.5 + aa[1] = 4. + aa[-2] = 4. + aa[-1] = 0.5 + + coef = compact.stdlu(aa) + + coef_verify = self.py_stdlu(aa) + + print ("RelError Norm: ", np.linalg.norm((coef - coef_verify)/coef_verify)) + + + def verify_nonp_lud2(self): + + print ("Non-periodic coef second derivative") + + nx = 512 + aa = np.ones(nx) * 3. + aa[0] = 2./11. + aa[1] = 10. + aa[-2] = 10. + aa[-1] = 2./11. + + coef = compact.stdlu(aa) + + coef_verify = self.py_stdlu(aa) + + print ("RelError Norm: ", np.linalg.norm((coef - coef_verify)/coef_verify)) + + + def py_stdlu(self, aa): + coef = np.ones(aa.shape)/aa[0] + + print ("coef.size = ", coef.size) + + for i in range(1,coef.size): + coef[i]=1.0/(aa[i]-coef[i-1]) + + return coef + + + def ddx (self, src): + + if src.shape != self.shape: + print ("error") + + nz, ny, nx = self.shape + + xsrc = np.zeros((ny, nx,), dtype=np.float64, order="F") + + # dst = np.zeros((nx, ny, nz,), order="F") + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.px: # Periodic BC + for i in range(nz): + dst[i] = compact.dfp(self.hx, src[i], 1) + + else: + for i in range(nz): + dst[i] = compact.dfnonp(self.hx, src[i], 1) + + # return np.swapaxes(dst, 1, 2) + return dst + + def ddy (self, src): + + if src.shape != self.shape: + print ("error") + + nz, ny, nx = self.shape + + #xsrc = np.zeros((ny, nx,), dtype=np.float64, order="F") + + # dst = np.zeros((nx, ny, nz,), order="F") + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.py: # Periodic BC + for i in range(nz): + dst[i] = compact.dfp(self.hx, src[i].T, 2).T + + else: + for i in range(nz): + dst[i] = compact.dfnonp(self.hx, src[i].T, 2).T + + # return np.swapaxes(dst, 1, 2) + return dst + + def ddz (self, src): + + if src.shape != self.shape: + print ("error") + + nz, ny, nx = self.shape + + + # dst = np.zeros((nx, ny, nz,), order="F") + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.pz: # Periodic BC + for i in range(ny): + dst[:,i,:] = compact.dfp(self.hx, src[:,i,:], 3) + + else: + for i in range(ny): + dst[:,i,:] = compact.dfnonp(self.hx, src[:,i,:], 3) + + # return np.swapaxes(dst, 1, 2) + return dst + + def port_nonp_coef (self): + + # SUBROUTINE nonp_lud(xyz,xx) + nz, ny, nx = self.shape + xx = nx + + lxf = np.zeros(xx) + lxs = np.zeros(xx) + + aa = np.zeros(xx) + aa[:] = 3. + + aa[0]=0.5 + aa[1]=4. + aa[-2]=4. + aa[-1]=0.5 + + # first derivative + compact.stdlu(aa,lxf) + + aa[:] = 5.5 + + aa[0]=2./11. + aa[1]=10. + aa[-2]=10. + aa[-1]=2./11. + + # second derivative + compact.stdlu(aa,lxs) + + compact.lxf = lxf + compact.lxs = lxs + + + + +def read_old_data (fname): + import struct + import sys + import os + with open(fname, 'rb') as f1 : + f1.seek(0) + + raw_info = f1.read(4+8*6+4)[4:-4] + t = struct.unpack('d', raw_info[ 0: 8])[0] + nx = struct.unpack('q', raw_info[ 8:16])[0] + ny = struct.unpack('q', raw_info[16:24])[0] + nz = struct.unpack('q', raw_info[24:32])[0] + count = nx*ny*nz + bSize = count*8 # size in bytes for a variable + + dummy_len = (4+8*3+4) + (4+8*2+4) + (4+8*2+4) + (4+8*2+4) + 4 + dummy = f1.read(dummy_len) + #dummy = f1.read(4) + + print (t, nx, ny, nz) + + #raw_field = f1.read(4+bSize*5+4)[4:-4] + V = np.fromfile(f1, dtype=np.float64, count=(3*count)).reshape((3,nz,ny,nx)) + s = np.fromfile(f1, dtype=np.float64, count=(2*count)).reshape((nz,ny,nx,2)) + + print (V.order) + print (s.order) + + print (V.shape) + print (s.shape) + + V.order="F" + s.order="F" + + print (V.shape) + print (s.shape) + + u = V[0] + v = V[1] + w = V[2] + + Y0 = s.T[0].T + Y1 = s.T[1].T + + return t, nx, ny, nz, u, v, w, Y0, Y1 + +def read_data (fname): + import struct + import sys + import os + with open(fname, 'rb') as f1 : + f1.seek(0) + + raw_info = f1.read(4+8*6+4)[4:-4] + t = struct.unpack('d', raw_info[ 0: 8])[0] + nx = struct.unpack('q', raw_info[ 8:16])[0] + ny = struct.unpack('q', raw_info[16:24])[0] + nz = struct.unpack('q', raw_info[24:32])[0] + count = nx*ny*nz + bSize = count*8 # size in bytes for a variable + + dummy_len = (4+8*3+4) + (4+8*2+4) + (4+8*2+4) + (4+8*2+4) + 4 + dummy = f1.read(dummy_len) + #dummy = f1.read(4) + + #raw_field = f1.read(4+bSize*5+4)[4:-4] + V = np.fromfile(f1, dtype=np.float64, count=(3*count)).reshape((3,nz,ny,nx)) + s = np.fromfile(f1, dtype=np.float64, count=(2*count)).reshape((2,nz,ny,nx)) + + print (V.flags) + print (s.flags) + + print (V.shape) + print (s.shape) + + u = V[0] + v = V[1] + w = V[2] + + Y0 = s[0] + Y1 = s[1] + + return t, nx, ny, nz, u, v, w, Y0, Y1 + + +def validate_trigonometric(): + + writeToFile = False + + shape = (256, 256, 512) + + nz, ny, nx = shape + + pi8 = np.arccos(-1.) + + l_0 = 2.0 + hyp=l_0*pi8/ny + hxp=hyp + hzp=hyp + + + cs = CompactScheme(nx, ny, nz, False, True, True, 4., 2., 2.) + + cs.test_ludcmp() + + cs.verify_nonp_lud1() + + cs.verify_nonp_lud2() + + cs.test_dfnonp() + + + + print ("Test ddx") + + Y1 = np.zeros(shape) + + true = np.zeros(shape) + + XX = np.arange(nx) * hxp + YY = np.arange(ny) * hyp + ZZ = np.arange(nz) * hzp + + + + print ("1-D sine test") + cos_fortran = compact.dfnonp(hxp, np.sin(1.1*XX).reshape((1,-1)), 1) + cos_exact = 1.1 * np.cos(1.1*XX).reshape((1,-1)) + + print ("Compact Scheme: ", cos_fortran.min(), cos_fortran.max()) + + print ("Exact : ",cos_exact.min(), cos_exact.max()) + + print ("Norm of relative errors: ", np.linalg.norm((cos_fortran - cos_exact)/cos_exact)) + # print (((cos_fortran - cos_exact)/cos_exact)) + + + + print ("3-D trigonometric test") + + zz, yy, xx = np.meshgrid(ZZ, YY, XX) + + Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] + + + def compare_3d_result (true1, dY1): + + print ("Calculated Min/Max", dY1.min(), dY1.max()) + print ("True Min/Max", true1.min(), true1.max()) + + eps = np.finfo(true1.dtype).eps + + relerr = (dY1 - true1) / (true1 + eps) + print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) + + print(" DDX Test ") + true[:] = (1.1 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + dY1 = cs.ddx(Y1)[:] + + compare_3d_result(true, dY1) + + + print(" DDY Test ") + Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] + true[:] = (-3.0 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + dY1 = cs.ddy(Y1)[:] + + compare_3d_result(true, dY1) + + + print(" DDZ Test ") + Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] + true[:] = (-2.0 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + dY1 = cs.ddz(Y1)[:] + + compare_3d_result(true, dY1) + + + if writeToFile: + y = np.memmap("phi", dtype=np.float64, mode="w+", shape=cs.shape) + y[:] = Y1[:] + + dydxtrue = np.memmap("dphitrue", dtype=np.float64, mode="w+", shape=cs.shape) + dydxtrue[:] = true1[:] + + dydx = np.memmap("dphi", dtype=np.float64, mode="w+", shape=cs.shape) + dydx[:] = dY1[:] + + + # cs.verify_nonp_coef() + + + +def test_dns_data(): + import sys + + file_name = "./fort.1000" + + answer = "./fort.2000" + + t, nx, ny, nz, u, v, w, Y0, Y1 = read_data(file_name) + + with open(answer, 'rb') as ans_file: + ans_file.seek(4) + ddx_answer = np.fromfile(ans_file, dtype=np.float64, count=Y1.size, ).reshape(Y1.shape) + + cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2) + + ddx = cs.ddx(Y1) + + print (ddx.min(), ddx.max()) + print (ddx_answer.min(), ddx_answer.max()) + relerr = (ddx - ddx_answer) / (ddx_answer) + err = (ddx - ddx_answer) + print ("Absolute Error", np.nanmin(err), np.nanmax(err)) + print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) + + + +if __name__ == "__main__": + + # validate_trigonometric() + + + print("DNS Field Test") + test_dns_data() +