research/incomp-flame-post
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

split ludcmp alloc and calc, explicit real type size

Browse source
This commit is contained in:
ignis 2019-07-06 02:43:46 +09:00
parent a71a854a0a
commit e282f210f6
1 changed files with 135 additions and 43 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

178
code/Compact.f90
View file

@ -1,14 +1,11 @@
MODULE Compact MODULE Compact
IMPLICIT NONE IMPLICIT NONE
PRIVATE REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, & lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs lzf,lzs,wzf,wzs
! lyzf,lyzs,wyzf,wyzs ! lyzf,lyzs,wyzf,wyzs
INTEGER :: nxc,nyc,nzc INTEGER :: nxc,nyc,nzc
REAL, PARAMETER :: ezero = 1.0e-14 REAL*8, PARAMETER :: ezero = 1.0e-14
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
CONTAINS CONTAINS
@ -20,6 +17,22 @@
nxc=nx nxc=nx
nyc=ny nyc=ny
nzc=nz 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' ! IF(nyc /= nzc) PRINT*,'ny should be equal nz'
! xp, yp, zp = 0 : periodic ! xp, yp, zp = 0 : periodic
@ -32,9 +45,6 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wxs(nxc),STAT=ierr) ALLOCATE(wxs(nxc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(1,nxc)
ELSE
CALL nonp_lud(1,nxc)
ENDIF ENDIF
ALLOCATE(lyf(nyc),STAT=ierr) ALLOCATE(lyf(nyc),STAT=ierr)
@ -46,9 +56,6 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wys(nyc),STAT=ierr) ALLOCATE(wys(nyc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
CALL p_lud(2,nyc)
ELSE
call nonp_lud(2,nyc)
ENDIF ENDIF
ALLOCATE(lzf(nzc),STAT=ierr) ALLOCATE(lzf(nzc),STAT=ierr)
@ -60,19 +67,95 @@
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
ALLOCATE(wzs(nzc),STAT=ierr) ALLOCATE(wzs(nzc),STAT=ierr)
IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' 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) CALL p_lud(3,nzc)
ELSE ELSE
call nonp_lud(3,nzc) call nonp_lud(3,nzc)
ENDIF ENDIF
! CALL x_lud END SUBROUTINE ludcmp_calculate
! CALL yz_lud
END SUBROUTINE ludcmp
SUBROUTINE nonp_lud(xyz,xx) SUBROUTINE nonp_lud(xyz,xx)
INTEGER :: i,xyz,xx INTEGER :: i,xyz,xx
REAL, DIMENSION(xx) :: aa REAL*8, DIMENSION(xx) :: aa
aa=3. aa=3.
aa(1)=0.5 ; aa(2)=4. aa(1)=0.5 ; aa(2)=4.
aa(xx-1)=4. ; aa(xx)=0.5 aa(xx-1)=4. ; aa(xx)=0.5
@ -92,7 +175,7 @@
SUBROUTINE p_lud(xyz,xx) SUBROUTINE p_lud(xyz,xx)
INTEGER :: i,xyz,xx INTEGER :: i,xyz,xx
REAL :: a REAL*8 :: a
a=3. ! first derivative a=3. ! first derivative
IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction 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.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
@ -105,8 +188,8 @@
SUBROUTINE stdlu(a,n,l) SUBROUTINE stdlu(a,n,l)
INTEGER :: n INTEGER :: n
REAL :: a(n),l(n) REAL*8 :: a(n),l(n)
REAL :: d REAL*8 :: d
INTEGER :: i INTEGER :: i
l(1)=1.0/a(1) l(1)=1.0/a(1)
DO i=2,n DO i=2,n
@ -117,9 +200,9 @@
SUBROUTINE ptdlu(a,n,l,w) SUBROUTINE ptdlu(a,n,l,w)
INTEGER :: n INTEGER :: n
REAL :: a,l(n),w(n) REAL*8 :: a,l(n),w(n)
INTEGER :: i INTEGER :: i
REAL :: aa(n),d REAL*8 :: aa(n),d
DO i=1,n-1 DO i=1,n-1
aa(i)=a aa(i)=a
@ -143,11 +226,14 @@
SUBROUTINE dfnonp(n,h,x,dx,nd,dir) SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h REAL*8,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 REAL*8 :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
! print *, "dfnonp received (nd,n)", nd, n
h1=1./h h1=1./h
@ -185,11 +271,14 @@
SUBROUTINE dfp(n,h,x,dx,nd,dir) SUBROUTINE dfp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h REAL*8,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL :: r1,r2,h1 REAL*8 :: r1,r2,h1
! print *, "dfnonp received (nd,n)", nd, n
h1=1./h h1=1./h
r1=7./3. r1=7./3.
@ -221,10 +310,10 @@
SUBROUTINE ptdslv(r,n,l,w,nd) SUBROUTINE ptdslv(r,n,l,w,nd)
INTEGER,INTENT(IN) :: n,nd INTEGER,INTENT(IN) :: n,nd
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l,w REAL*8,INTENT(IN),DIMENSION(:) :: l,w
INTEGER i,j INTEGER i,j
REAL, DIMENSION(nd) :: sum REAL*8, DIMENSION(nd) :: sum
DO j=1,nd DO j=1,nd
sum(j)=w(1)*r(j,1) sum(j)=w(1)*r(j,1)
r(j,1)=r(j,1)*l(1) r(j,1)=r(j,1)*l(1)
@ -249,11 +338,13 @@
SUBROUTINE d2fp(n,h,x,dx,nd,dir) SUBROUTINE d2fp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h REAL*8,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL :: h2,r1,r2,t1,t2 REAL*8 :: h2,r1,r2,t1,t2
h2=1./(h*h) h2=1./(h*h)
r1=6. r1=6.
r2=3./8. r2=3./8.
@ -312,10 +403,10 @@
SUBROUTINE tdslv(r,n,l,nd) SUBROUTINE tdslv(r,n,l,nd)
INTEGER,INTENT(IN) :: n,nd INTEGER,INTENT(IN) :: n,nd
REAL,INTENT(INOUT),DIMENSION(nd,n) :: r REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l REAL*8,INTENT(IN),DIMENSION(:) :: l
INTEGER i,j INTEGER i,j
REAL t1 REAL*8 t1
DO j=1,nd DO j=1,nd
r(j,1)=r(j,1)*l(1) r(j,1)=r(j,1)*l(1)
ENDDO ENDDO
@ -334,11 +425,12 @@
SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h REAL*8,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j INTEGER :: i,j
REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2 REAL*8 :: h2,r1,r2,r3,a,b,c,e,t1,t2
h2=1./(h*h) h2=1./(h*h)
r1=6. r1=6.