Permalink
Browse files

Faster implementation of implicit routines.

  • Loading branch information...
stevensrjam committed Feb 17, 2016
1 parent 042737a commit 2f6ab93e2f11d57cce6326faf2d8eb00ac4c3cf4
Showing with 52 additions and 49 deletions.
  1. +18 −23 SolveImpEqnUpdate_Temp.F90
  2. +18 −15 SolveImpEqnUpdate_X.F90
  3. +16 −11 SolveImpEqnUpdate_YZ.F90
View
@@ -13,8 +13,8 @@ subroutine SolveImpEqnUpdate_Temp
use local_arrays, only : temp,rhs
use decomp_2d, only: xstart,xend
implicit none
real, dimension(nx) :: amkl,apkl,ackl, fkl
integer :: jc,kc,info,ipkv(nx),ic
real, dimension(nx) :: amkl,apkl,ackl
integer :: jc,kc,info,ipkv(nx),ic,nrhs
real :: betadx,ackl_b
real :: amkT(nx-1),ackT(nx),apkT(nx-1),appk(nx-2)
@@ -45,30 +45,25 @@ subroutine SolveImpEqnUpdate_Temp
! No solving is done in this call.
call dgttrf(nx,amkT,ackT,apkT,appk,ipkv,info)
nrhs=(xend(3)-xstart(3)+1)*(xend(2)-xstart(2)+1)
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
! Normalize RHS of equation
fkl(1)= 0.d0
do kc=2,nxm
ackl_b=1.0/(1.-ac3ssk(kc)*betadx)
fkl(kc)=rhs(kc,jc,ic)*ackl_b
end do
fkl(nx)= 0.d0
! Solve equation using LAPACK library
call dgttrs('N',nx,1,amkT,ackT,apkT,appk,ipkv,fkl,nx,info)
! Update temperature field
do kc=2,nxm
temp(kc,jc,ic) = temp(kc,jc,ic) + fkl(kc)
do jc=xstart(2),xend(2)
do kc=2,nxm
ackl_b=1.0/(1.0-ac3ssk(kc)*betadx)
rhs(kc,jc,ic)=rhs(kc,jc,ic)*ackl_b
end do
end do
end do
call dgttrs('N',nx,nrhs,amkT,ackT,apkT,appk,ipkv,rhs,nx,info)
enddo
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
do kc=2,nxm
temp(kc,jc,ic)=temp(kc,jc,ic) + rhs(kc,jc,ic)
end do
end do
end do
return
View
@@ -14,11 +14,11 @@ subroutine SolveImpEqnUpdate_X
use local_arrays, only : vx,rhs
use decomp_2d, only: xstart,xend
implicit none
real, dimension(nx) :: amkl,apkl,ackl, fkl
real, dimension(nx) :: amkl,apkl,ackl
real :: amkT(nx-1),apkT(nx-1)
real :: appk(nx-2)
real :: ackT(nx)
integer :: jc,kc,info,ic
integer :: jc,kc,info,ic,nrhs
integer :: ipkv(nx)
real :: betadx,ackl_b
@@ -43,21 +43,24 @@ subroutine SolveImpEqnUpdate_X
call dgttrf(nx,amkT,ackT,apkT,appk,ipkv,info)
nrhs=(xend(3)-xstart(3)+1)*(xend(2)-xstart(2)+1)
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
fkl(1)= 0.d0
do kc=2,nxm
ackl_b=1.0d0/(1.0d0-ac3ck(kc)*betadx)
fkl(kc) = rhs(kc,jc,ic)*ackl_b
enddo
fkl(nx)= 0.d0
call dgttrs('N',nx,1,amkT,ackT,apkT,appk,ipkv,fkl,nx,info)
do jc=xstart(2),xend(2)
do kc=2,nxm
ackl_b=1.0/(1.0-ac3ck(kc)*betadx)
rhs(kc,jc,ic)=rhs(kc,jc,ic)*ackl_b
end do
end do
end do
call dgttrs('N',nx,nrhs,amkT,ackT,apkT,appk,ipkv,rhs,nx,info)
do kc=2,nxm
vx(kc,jc,ic)=vx(kc,jc,ic) + fkl(kc)
enddo
enddo
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
do kc=2,nxm
vx(kc,jc,ic)=vx(kc,jc,ic) + rhs(kc,jc,ic)
end do
end do
end do
return
View
@@ -17,8 +17,8 @@ subroutine SolveImpEqnUpdate_YZ(q,rhs)
& xstart(3)-lvlhalo:xend(3)+lvlhalo)
real,intent(inout) :: rhs(1:nx,xstart(2):xend(2), &
& xstart(3):xend(3))
real, dimension(nx) :: amkl,apkl,ackl,fkl
integer :: jc,kc,info,ipkv(nxm),ic
real, dimension(nx) :: amkl,apkl,ackl
integer :: jc,kc,info,ipkv(nxm),ic,nrhs
real :: betadx,ackl_b
real :: amkT(nxm-1),ackT(nxm),apkT(nxm-1),appk(nx-3)
@@ -37,19 +37,24 @@ subroutine SolveImpEqnUpdate_YZ(q,rhs)
call dgttrf(nxm,amkT,ackT,apkT,appk,ipkv,info)
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
do kc=1,nxm
ackl_b=1.0d0/(1.-ac3sk(kc)*betadx)
fkl(kc)=rhs(kc,jc,ic)*ackl_b
nrhs=(xend(3)-xstart(3)+1)*(xend(2)-xstart(2)+1)
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
do kc=1,nxm
ackl_b=1.0/(1.0-ac3sk(kc)*betadx)
rhs(kc,jc,ic)=rhs(kc,jc,ic)*ackl_b
end do
end do
end do
call dgttrs('N',nxm,1,amkT,ackT,apkT,appk,ipkv,fkl,nxm,info)
call dgttrs('N',nxm,nrhs,amkT,ackT,apkT,appk,ipkv,rhs,nx,info)
do kc=1,nxm
q(kc,jc,ic)=q(kc,jc,ic) + fkl(kc)
do ic=xstart(3),xend(3)
do jc=xstart(2),xend(2)
do kc=1,nxm
q(kc,jc,ic)=q(kc,jc,ic) + rhs(kc,jc,ic)
end do
end do
end do
end do
return

0 comments on commit 2f6ab93

Please sign in to comment.