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

Initial commit

Browse source
This commit is contained in:
ignis 2019-04-01 15:12:36 +09:00
commit 7b4dec31f5
4 changed files with 2195 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

386
code/Compact.f90 Normal file
View file

@ -0,0 +1,386 @@
MODULE Compact
IMPLICIT NONE
PRIVATE
REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
lyf,lys,wyf,wys, &
lzf,lzs,wzf,wzs
! lyzf,lyzs,wyzf,wyzs
INTEGER :: nxc,nyc,nzc
REAL, PARAMETER :: ezero = 1.0e-14
PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
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
! 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'
CALL p_lud(1,nxc)
ELSE
CALL nonp_lud(1,nxc)
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'
CALL p_lud(2,nyc)
ELSE
call nonp_lud(2,nyc)
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'
CALL p_lud(3,nzc)
ELSE
call nonp_lud(3,nzc)
ENDIF
! CALL x_lud
! CALL yz_lud
END SUBROUTINE ludcmp
SUBROUTINE nonp_lud(xyz,xx)
INTEGER :: i,xyz,xx
REAL, 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 :: 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 :: a(n),l(n)
REAL :: d
INTEGER :: i
l(1)=1.0/a(1)
DO i=2,n
d=a(i)-l(i-1)
l(i)=1.0/d
ENDDO
END SUBROUTINE stdlu
SUBROUTINE ptdlu(a,n,l,w)
INTEGER :: n
REAL :: a,l(n),w(n)
INTEGER :: i
REAL :: 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 dfnonp(n,h,x,dx,nd,dir)
INTEGER,INTENT(IN) :: n,nd,dir
REAL,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4
h1=1./h
r1=7./3.
r2=1./12.
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
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,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: r1,r2,h1
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,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l,w
INTEGER i,j
REAL, 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,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: 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,INTENT(INOUT),DIMENSION(nd,n) :: r
REAL,INTENT(IN),DIMENSION(:) :: l
INTEGER i,j
REAL 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,INTENT(IN) :: h
REAL,INTENT(IN),DIMENSION(nd,n) :: x
REAL,INTENT(OUT),DIMENSION(nd,n) :: dx
INTEGER :: i,j
REAL :: 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

20
code/makefile Normal file
View file

@ -0,0 +1,20 @@
#flags = -r8 -i8 -msse3 -O3 -axT -m64 -fPIC -i_dynamic -openmp
#Neuman flags
#flags = -O3 -r8 -i8 -openmp -msse3 -axSSSE3 -m64 -fPIC -i_dynamic
#Cluster_2(16.161) flags
flags = -O3 -r8 -i8 -qopenmp -msse3 -axSSSE3 -m64 -fPIC
compiler = ifort
Post_Cbar_FSD_Incomp : test.o post.mod Compact.mod
${compiler} -o x-edge-cold-bc-uPrime-hybrid ${flags} test.o post.o Compact.o
test.o : test.f90 post.mod
${compiler} -c ${flags} test.f90
post.mod: post.f90 Compact.mod
${compiler} -c ${flags} post.f90
Compact.mod: Compact.f90
${compiler} -c ${flags} Compact.f90
clean:
rm *.o *.mod x-edge-cold-bc-uPrime-hybrid

1778
code/post.f90 Normal file
View file

File diff suppressed because it is too large Load diff

11
code/test.f90 Normal file
View file

@ -0,0 +1,11 @@
! This code has been written since 2011.01.08
! This code is written by Dongkyu Lee (THDHRKDG)
! Purpose of this code is to postprocess results of stagnating DNS code
PROGRAM test
USE post
USE Compact
CALL main
END PROGRAM