corse grained and term-wise derivative computation
This commit is contained in:
ignis
parent
6933297130
commit
749e77de32
1 changed files with 227 additions and 66 deletions
293
m_fdm_calc.f90
293
m_fdm_calc.f90
|
|
@ -10,8 +10,11 @@ module m_fdm_calc
|
|||
|
||||
real*8, dimension(:,:,:), allocatable :: u_,v_,w_
|
||||
real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
|
||||
real*8, dimension(:,:), allocatable :: fz, dfz, fzz, dfzz
|
||||
real*8, dimension(:,:), allocatable :: fzl, fzu, fzzl, fzzu
|
||||
real*8, dimension(:,:), allocatable :: fzzl, fzzu
|
||||
|
||||
real*8, dimension(:,:), allocatable :: fbuf1, fbuf2, fbuf3, fbuf4
|
||||
real*8, dimension(:,:), allocatable :: yxbuf1, yxbuf2, yxbuf3, yxbuf4
|
||||
real*8, dimension(:,:), allocatable :: xybuf1, xybuf2, xybuf3, xybuf4
|
||||
|
||||
real*8 :: in_yr,out_yr,refwr,minf
|
||||
integer :: fullsavenum !,svfx,svfy
|
||||
|
|
@ -240,17 +243,24 @@ module m_fdm_calc
|
|||
allocate(y2(nx,ny,nz,neq))
|
||||
allocate(yf(nx,ny,nz,neq))
|
||||
|
||||
allocate(fz(4*nx*ny,nz))
|
||||
allocate(dfz(4*nx*ny,nz))
|
||||
allocate(fzz(nx*ny,nz))
|
||||
allocate(dfzz(nx*ny,nz))
|
||||
|
||||
allocate(fzu(4*nx*ny,2))
|
||||
allocate(fzl(4*nx*ny,2))
|
||||
|
||||
allocate(fzzu(nx*ny,2))
|
||||
allocate(fzzl(nx*ny,2))
|
||||
|
||||
allocate(fbuf1(nx*ny,nz))
|
||||
allocate(fbuf2(nx*ny,nz))
|
||||
allocate(fbuf3(nx*ny,nz))
|
||||
allocate(fbuf4(nx*ny,nz))
|
||||
|
||||
allocate(xybuf1(nx,ny))
|
||||
allocate(xybuf2(nx,ny))
|
||||
allocate(xybuf3(nx,ny))
|
||||
allocate(xybuf4(nx,ny))
|
||||
|
||||
allocate(yxbuf1(ny,nx))
|
||||
allocate(yxbuf2(ny,nx))
|
||||
allocate(yxbuf3(ny,nx))
|
||||
allocate(yxbuf4(ny,nx))
|
||||
|
||||
y1=0.0
|
||||
y2=0.0
|
||||
yf=0.0
|
||||
|
|
@ -467,21 +477,18 @@ module m_fdm_calc
|
|||
|
||||
implicit none
|
||||
|
||||
integer :: i,j,k,xx,yy,zz,ii
|
||||
integer, parameter :: nb = 4
|
||||
|
||||
integer :: i,j,k,xx,yy,zz,ii,jj,kk
|
||||
integer :: n
|
||||
integer :: idx1, idx2
|
||||
real*8 :: wrate,yr,yp
|
||||
real*8 :: r1_(xx,yy,zz,neq),f_(xx,yy,zz,neq)
|
||||
real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
|
||||
real*8 :: ux(4,xx),dux(4,xx),d2ux(xx)
|
||||
real*8 :: uy(4,yy),duy(4,yy),d2uy(yy)
|
||||
real*8 :: uz(4,zz),duz(4,zz),d2uz(zz)
|
||||
real*8 :: uux(xx)
|
||||
real*8 :: uuy(yy)
|
||||
real*8 :: uuz(zz)
|
||||
real*8 :: y
|
||||
|
||||
|
||||
! reaction source term
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
|
|
@ -534,96 +541,147 @@ module m_fdm_calc
|
|||
|
||||
else
|
||||
|
||||
! -( d(rho*w)/dz )
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
idx2 = xx*(j-1)+i
|
||||
idx1 = (idx2-1)*4
|
||||
fz(idx1+1,k) = r1_(i,j,k,1) ! 1:rho
|
||||
fz(idx1+2,k) = r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y
|
||||
fz(idx1+3,k) = r1_(i,j,k,1)*ww_(i,j,k) ! 3:rho*w
|
||||
fz(idx1+4,k) = r1_(i,j,k,2)*ww_(i,j,k) ! 4:rho*w*Y
|
||||
fzz(idx2,k) = r1_(i,j,k,2)/r1_(i,j,k,1)
|
||||
!fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k) ! :rho*w*YO
|
||||
!fzz(idx2,k) = r1_(2,i,j,k)/r1_(1,i,j,k) ! :YO
|
||||
fbuf1(i+xx*(j-1),k) = r1_(i,j,k,1)*ww_(i,j,k) ! rho*w
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
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)
|
||||
CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3)
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
f_(i,j,k,1) = f_(i,j,k,1) - fbuf2(i+xx*(j-1),k) ! continuity
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
|
||||
! = -( d(rho*w*Yr)/dz )
|
||||
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
fbuf1(i+xx*(j-1),k) = r1_(i,j,k,2)*ww_(i,j,k) ! rho*w*Y
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3)
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
idx2 = xx*(j-1)+i
|
||||
idx1 = (idx2-1)*4
|
||||
! -( d(rho*w)/dz )
|
||||
f_(i,j,k,1) = f_(i,j,k,1) - dfz(idx1+3,k) ! continuity
|
||||
f_(i,j,k,2) = f_(i,j,k,2) - fbuf2(i+xx*(j-1),k) ! species conserv.
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
|
||||
! = -( d(rho*w*Yr)/dz )
|
||||
! + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
|
||||
f_(i,j,k,2) = f_(i,j,k,2) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2,k)+dfz(idx1+1,k)*dfz(idx1+2,k)) ! species conserv.
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
fbuf1(i+xx*(j-1),k) = r1_(i,j,k,2)/r1_(i,j,k,1) ! Y
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
CALL pdfp (fbuf1, fzzl, fzzu, fbuf2, hy, xx*yy, zz, yy, 3)
|
||||
CALL pd2fp(fbuf1, fzzl, fzzu, fbuf3, hy, xx*yy, zz, yy, 3)
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
fbuf1(i+xx*(j-1),k) = r1_(i,j,k,1) ! rho
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
CALL pdfp (fbuf1, fzzl, fzzu, fbuf4, hy, xx*yy, zz, yy, 3)
|
||||
|
||||
fbuf2 = fbuf2 * fbuf4
|
||||
fbuf1 = fbuf1 * fbuf3 + fbuf2
|
||||
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
f_(i,j,k,2) = f_(i,j,k,2) + diff*fbuf1(i+xx*(j-1),k) ! species conserv.
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
endif
|
||||
|
||||
|
||||
! x-direction
|
||||
DO k=1,zz
|
||||
DO j=1,yy
|
||||
DO i=1,xx
|
||||
ux(1,i)=r1_(i,j,k,1) ! 1:rho
|
||||
ux(2,i)=r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y
|
||||
ux(3,i)=ux(1,i)*uu_(i,j,k) ! 3:rho*u
|
||||
ux(4,i)=ux(3,i)*ux(2,i) ! 4:rho*u*Y
|
||||
uux (i)=ux(2,i)
|
||||
ENDDO
|
||||
|
||||
CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
|
||||
CALL d2fnonp(xx,hx,uux(:),d2ux(:),1,1)
|
||||
|
||||
DO i=1,xx
|
||||
|
||||
! -( d(rho*u)/dx )
|
||||
f_(i,j,k,1) = f_(i,j,k,1) - dux(3,i) ! continuity
|
||||
CALL tp2mul (r1_(:,:,k,1), uu_(:,:,k), yxbuf1, xx, yy) ! rho*u
|
||||
CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u)
|
||||
CALL tp2sub (f_(:,:,k,1), yxbuf2, xx, yy) ! continuity
|
||||
|
||||
! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
|
||||
! = -( d(rho*u*Yr)/dx )
|
||||
! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
|
||||
f_(i,j,k,2) = f_(i,j,k,2) - dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) ! species conservation
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
CALL tp2mul (r1_(:,:,k,2), uu_(:,:,k), yxbuf1, xx, yy) ! rho*u*Y
|
||||
CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u*Y)
|
||||
CALL tp2sub (f_(:,:,k,2), yxbuf2, xx, yy) ! species conservation
|
||||
|
||||
CALL tp2div (r1_(:,:,k,2), r1_(:,:,k,1), yxbuf1, xx, yy)! Y
|
||||
CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(Y)
|
||||
CALL d2fnonp(xx,hx,yxbuf1,yxbuf3,yy,1) ! d2/dx2(Y)
|
||||
|
||||
CALL tp2 (r1_(:,:,k,1), yxbuf1, xx, yy) ! rho
|
||||
CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho)
|
||||
|
||||
yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho)
|
||||
yxbuf1 = -diff*(yxbuf2 + yxbuf1 * yxbuf3) ! -D( ... + rho * d2/dx2(Y))
|
||||
|
||||
CALL tp2sub (f_(:,:,k,2), yxbuf1, xx, yy) ! species conservation
|
||||
|
||||
ENDDO
|
||||
|
||||
|
||||
!! y-direction
|
||||
DO k=1,zz
|
||||
DO i=1,xx
|
||||
DO j=1,yy
|
||||
uy(1,j)=r1_(i,j,k,1) ! 1:rho
|
||||
uy(2,j)=r1_(i,j,k,2)/r1_(i,j,k,1) ! 2:Y
|
||||
uy(3,j)=uy(1,j)*vv_(i,j,k) ! 3:rho*v
|
||||
uy(4,j)=uy(3,j)*uy(2,j) ! 4:rho*v*Y
|
||||
uuy (j)=uy(2,j)
|
||||
ENDDO
|
||||
|
||||
CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
|
||||
CALL d2fp(yy,hy,uuy(:),d2uy(:),1,2)
|
||||
! -( d(rho*v)/dy )
|
||||
|
||||
DO j=1,yy
|
||||
! -( d(rho*v)/dy )
|
||||
f_(i,j,k,1)=f_(i,j,k,1)-duy(3,j) ! continuity
|
||||
xybuf1(:,:)=r1_(:,:,k,1)*vv_(:,:,k) ! rho*v
|
||||
CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(rho*v)
|
||||
|
||||
f_(:,:,k,1)=f_(:,:,k,1)-xybuf2(:,:) ! continuity
|
||||
|
||||
! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
|
||||
! = -( d(rho*v*Yr)/dy )
|
||||
! + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
|
||||
f_(i,j,k,2)=f_(i,j,k,2)-duy(4,j) + diff*(uy(1,j)*d2uy(j)+duy(1,j)*duy(2,j)) ! species conserv.
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
xybuf1(:,:) = r1_(:,:,k,2)*vv_(:,:,k) ! rho*v*Y
|
||||
CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(rho*v*Y)
|
||||
|
||||
f_(:,:,k,2) = f_(:,:,k,2) - xybuf2(:,:) ! species conserv.
|
||||
|
||||
xybuf1(:,:) = r1_(:,:,k,2)/r1_(:,:,k,1) ! Y
|
||||
CALL dfp(yy,hy,xybuf1,xybuf2,xx,2) ! d/dy(Y)
|
||||
CALL d2fp(yy,hy,xybuf1,xybuf3,xx,2) ! d2/dy2(Y)
|
||||
|
||||
xybuf1(:,:) = r1_(:,:,k,1) ! rho
|
||||
CALL dfp(yy,hy,xybuf1,xybuf4,xx,2) ! d/dy(rho)
|
||||
|
||||
xybuf2 = xybuf2 * xybuf4 ! d/dy(Y) * d/dy(rho)
|
||||
xybuf1 = xybuf2 + xybuf1 * xybuf3 ! ... + rho * d2/dy2(Y)
|
||||
|
||||
f_(:,:,k,2) = f_(:,:,k,2) + diff*xybuf1(:,:) ! species conserv.
|
||||
|
||||
ENDDO
|
||||
|
||||
|
||||
|
|
@ -735,6 +793,109 @@ module m_fdm_calc
|
|||
return
|
||||
END SUBROUTINE substep
|
||||
|
||||
subroutine tp2 (a, b, n1, n2)
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: n1, n2
|
||||
real*8,intent(in) :: a(n1,n2)
|
||||
real*8,intent(out) :: b(n2,n1)
|
||||
integer :: i,j,ii,jj
|
||||
integer,parameter :: nb = 16
|
||||
|
||||
DO jj=1,n2,nb
|
||||
DO ii=1,n1,nb
|
||||
|
||||
DO j=jj,jj+nb-1
|
||||
DO i=ii,ii+nb-1
|
||||
b(j,i) = a(i,j)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
end subroutine tp2
|
||||
|
||||
|
||||
subroutine tp2div (a, b, c, n1, n2)
|
||||
! c = a / b
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: n1, n2
|
||||
real*8,intent(in) :: a(n1,n2), b(n1,n2)
|
||||
real*8,intent(out) :: c(n2,n1)
|
||||
integer :: i,j,ii,jj
|
||||
integer,parameter :: nb = 16
|
||||
|
||||
DO jj=1,n2,nb
|
||||
DO ii=1,n1,nb
|
||||
|
||||
DO j=jj,jj+nb-1
|
||||
DO i=ii,ii+nb-1
|
||||
c(j,i) = a(i,j) / b(i,j)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
|
||||
end subroutine tp2div
|
||||
|
||||
|
||||
subroutine tp2mul (a, b, c, n1, n2)
|
||||
! c = a * b
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: n1, n2
|
||||
real*8,intent(in) :: a(n1,n2), b(n1,n2)
|
||||
real*8,intent(out) :: c(n2,n1)
|
||||
integer :: i,j,ii,jj
|
||||
integer,parameter :: nb = 16
|
||||
|
||||
DO jj=1,n2,nb
|
||||
DO ii=1,n1,nb
|
||||
|
||||
DO j=jj,jj+nb-1
|
||||
DO i=ii,ii+nb-1
|
||||
c(j,i) = a(i,j) * b(i,j)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
end subroutine tp2mul
|
||||
|
||||
|
||||
subroutine tp2sub (a, b, n1, n2)
|
||||
! a = a - transpose(b)
|
||||
|
||||
implicit none
|
||||
|
||||
integer,intent(in) :: n1, n2
|
||||
real*8,intent(inout) :: a(n1,n2)
|
||||
real*8,intent(in) :: b(n2,n1)
|
||||
integer :: i,j,ii,jj
|
||||
integer,parameter :: nb = 16
|
||||
|
||||
DO jj=1,n2,nb
|
||||
DO ii=1,n1,nb
|
||||
|
||||
DO j=jj,jj+nb-1
|
||||
DO i=ii,ii+nb-1
|
||||
a(i,j) = a(i,j) - b(j,i)
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
ENDDO
|
||||
ENDDO
|
||||
|
||||
end subroutine tp2sub
|
||||
|
||||
|
||||
end module m_fdm_calc
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue