diff options
Diffstat (limited to 'src/Stab3d.F')
| -rw-r--r-- | src/Stab3d.F | 280 |
1 files changed, 280 insertions, 0 deletions
diff --git a/src/Stab3d.F b/src/Stab3d.F new file mode 100644 index 0000000..1aaa7da --- /dev/null +++ b/src/Stab3d.F @@ -0,0 +1,280 @@ + subroutine bicgst3d(cc,ce,cw,cn,cs,ct,cb,u,rhs, + & eps,rmax,ier,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c + implicit none +c + integer,intent(in) :: im,jm,km + real*8,intent(inout) :: cc(im,jm,km),cn(im,jm,km),cs(im,jm,km), + $ ce(im,jm,km),cw(im,jm,km),ct(im,jm,km),cb(im,jm,km) + real*8,intent(out) :: eps + real*8,intent(out) :: rmax + real*8 :: u(im,jm,km),rhs(im,jm,km) +c Local variable + integer ncyc + integer iscale,i,j,k,ier +c +c****************************************************************************** +c + + ncyc = (im-2)*(jm-2)*(km-2) + + ier=0 + +* +* Determine whether we can diagonally scale the problem to speed +* convergence. Can only be done if there are no zeros on the main +* diagonal (ie. central difference coefficient). +* + iscale=1 + do k = 2,km-1 + do j = 2,jm-1 + do i = 2,im-1 + if (cc(i,j,k).eq.0.) iscale = 0 + enddo + enddo + enddo + +* +* Do the diagonal scaling if we can. +* + if (iscale.eq.1) then + + + do k = 2,km-1 + do j = 2,jm-1 + do i = 2,im-1 + rhs(i,j,k)=rhs(i,j,k)/cc(i,j,k) + cb(i,j,k)=cb(i,j,k)/cc(i,j,k) + ct(i,j,k)=ct(i,j,k)/cc(i,j,k) + cw(i,j,k)=cw(i,j,k)/cc(i,j,k) + ce(i,j,k)=ce(i,j,k)/cc(i,j,k) + cs(i,j,k)=cs(i,j,k)/cc(i,j,k) + cn(i,j,k)=cn(i,j,k)/cc(i,j,k) + cc(i,j,k)=1. + enddo + enddo + enddo + + endif + +* +* Now call the bicgstab routine +* + + call bicgstab (cc,cn,cs,ce,cw,ct,cb,u,rhs,eps,ncyc,rmax,ier, + & im,jm,km) + if (rmax.gt.eps) then + ier = -1 + print*,'Did not converge' + print*,' maximum residual = ',rmax + print*,' tolerance = ',eps + endif + + return + end +c +c****************************************************************** +c + subroutine bicgstab (cc,cn,cs,ce,cw,ct,cb,x,r,tol,ncyc,rnorm, + $ ier,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c + implicit none +c + integer,intent(in) :: ncyc,im,jm,km + real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km) + real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km) + real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km) + real*8,intent(out) :: cb(im*jm*km) + real*8,intent(out) :: tol,rnorm + integer,intent(out) :: ier + real*8 x(im*jm*km),r(im*jm*km) +c Local variables + integer :: i,j,k,kk + real*8 :: p(im*jm*km),Ap(im*jm*km),w(im*jm*km),As(im*jm*km) + real*8 :: omega, chi,chi1,chi2, delta, deltap, pp +* +*********************************************************************** +* + + do i = 1,im*jm*km + p(i) = 0. + Ap(i) = 0. + w(i) = 0. + As(i) = 0. + enddo + + kk = 0 + 10 call usermv(cc,cn,cs,ce,cw,ct,cb,x,Ap,im,jm,km) + + + do i = 1,im*jm*km + r(i) = r(i)-Ap(i) + p(i) = r(i) + enddo + +c delta = sum(r) + delta = 0. + + do i = 1,im*jm*km + delta = delta+r(i) + enddo + + if (delta.eq.0.) then + ier=-1 + return + endif + + call usermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km) + +c phi = sum(Ap) + pp = 0. + + do i = 1,im*jm*km + pp = pp+Ap(i) + enddo + pp = pp/delta + + if (pp.eq.0.) then + ier=-1 + return + endif + +c rnorm = sum(r**2) + rnorm = 0. + + do i = 1,im*jm*km + rnorm = rnorm + r(i)**2 + enddo + rnorm=sqrt(rnorm) + +c Test if initial guess is great (residual less than tolerance) + if (rnorm .lt. tol) return + + 1 continue + kk = kk + 1 + + omega = 1./pp + + + do i = 1,im*jm*km + w(i) = r(i) - omega*Ap(i) + enddo + + call usermv(cc,cn,cs,ce,cw,ct,cb,w,As,im,jm,km) + +c chi1 = sum(As*w) + chi1 = 0. + + do i = 1,im*jm*km + chi1 = chi1+As(i)*w(i) + enddo +c chi2 = sum(As**2) + chi2 = 0. + + do i = 1,im*jm*km + chi2 = chi2+As(i)**2 + enddo + + chi = chi1/chi2 + + + do i = 1,im*jm*km + r(i) = w(i) - chi*As(i) + x(i) = x(i) + omega*p(i) + chi*w(i) + enddo + + deltap = delta +c delta = sum(r) + delta = 0. + + do i = 1,im*jm*km + delta = delta+r(i) + enddo + + if (delta.eq.0.) then + goto 10 + endif + + + do i = 1,im*jm*km + p(i) = r(i) + (p(i)-chi*Ap(i))*omega* + & delta/(deltap*chi) + enddo + + call usermv(cc,cn,cs,ce,cw,ct,cb,p,Ap,im,jm,km) + +c phi = sum(Ap) + pp = 0. + + do i = 1,im*jm*km + pp = pp+Ap(i) + enddo + pp=pp/delta + + if (pp.eq.0.) then + goto 10 + endif + + if (kk.gt.ncyc) then + print*,' BI-CGStab solver reached maximum nuber of iterations.' + ier=-1 + return + endif + +c rnorm = sum(r**2) + rnorm = 0. + + do i = 1,im*jm*km + rnorm = rnorm+r(i)**2 + enddo + rnorm=sqrt(rnorm) + + if (rnorm .gt. tol) goto 1 + print*,'-----------' + write(*,*) "ncyc=",kk + print*,'-----------------------------------------------------------' + return + end +c +c****************************************************************** +c + subroutine usermv(cc,cn,cs,ce,cw,ct,cb,x,y,im,jm,km) +c +c This routine was lifted from stab.f. Minor modifications have +c been made. +c +c Be careful that the c's are zero on their outer boundary!! +c + implicit none +c + integer,intent(in) :: im,jm,km + real*8,intent(out) :: cc(im*jm*km),cn(im*jm*km) + real*8,intent(out) :: cs(im*jm*km),ce(im*jm*km) + real*8,intent(out) :: cw(im*jm*km),ct(im*jm*km) + real*8,intent(out) :: cb(im*jm*km) + real*8 :: x(im*jm*km), y(im*jm*km) +c Local variables + integer :: i, j, k +* +*********************************************************************** +* +* + do i = im*jm+im+1,im*(jm*km-jm-1) + y(i) = cw(i)*x(i-1) + & +cc(i)*x(i) + & +ce(i)*x(i+1) + & +cn(i)*x(i+im) + & +cs(i)*x(i-im) + & +ct(i)*x(i+im*jm) + & +cb(i)*x(i-im*jm) + enddo +c + return + end + |