diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 0cce4b4..6d63cfd 100644 --- a/m_fdm_calc.f90 +++ b/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