added extract ddxc, timestamp to log
numpy memmap py3 compatibility fix
This commit is contained in:
ignis
parent
dba0c40891
commit
d947ffe62a
6 changed files with 1120 additions and 12 deletions
539
Compact.f90
Normal file
539
Compact.f90
Normal file
|
|
@ -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
|
||||
2
compact_doc.py
Normal file
2
compact_doc.py
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
from compact import compact
|
||||
print compact.__doc__
|
||||
15
extract_c.py
15
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()
|
||||
|
|
|
|||
33
extract_ddxc.py
Normal file
33
extract_ddxc.py
Normal file
|
|
@ -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()
|
||||
15
extract_u.py
15
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()
|
||||
|
|
|
|||
528
pycompact.py
Normal file
528
pycompact.py
Normal file
|
|
@ -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()
|
||||
|
||||
Loading…
Add table
Reference in a new issue