[Nek5000-users] Residuals for post-processing
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Wed Apr 21 09:27:12 CDT 2010
Hi Markus,
pressure residuals are in gmres
You could output them via:
call outpost(x,y,z,res,t,'res')
where x,y,z,t are arrays of your choice and res
is the array you want to look at.
Nek will convert the res mesh to the velocity mesh
if needed. I usually follow such a call with
call exitti('quit in gmres$',n)
because prepost() uses some scratch common blocks
that are potentially shared with the calling routine
(not in the case of gmres, I think...) -- so data
might be corrupted after such a call. The idea of
prepost was to be called in between timesteps, when
such memory re-use would be admissable. The alternative
would be to define your own storage in a common block,
save the residuals to this locale, then loop through
calls to outpost in userchk afterwards... I've done
that in the past.
All this being said, I'm not 100% certain that the
iterative solvers is the place to look for the difficulty.
Most problems arise from spatial resolution issues...or
temporal stability problems. The latter can be identified
by drastic reduction in dt.
What exactly is blowing up for you?
Paul
On Wed, 21 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:
> Hi,
>
> I'd like to save the residuals for the pressure iterations (PN/PN-2 as well
> as PN/PN) into a scalar variable (scalar 0, instead of temperature say) for
> each grid point in order to see where my solution diverges.
> Is there a way to add something to the .usr file to do that? I couldn't find
> the variable in the code and also wouldn't know how to save something from
> the lx2=lx1-2 grid onto the lx1 grid.
> Would this also be possible with the velocity residuals?
>
> Thanks,
> Markus
>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>
uuuu
c-----------------------------------------------------------------------
subroutine uzawa_gmres(res,h1,h2,h2inv,intype,iter)
c Solve the pressure equation by right-preconditioned
c GMRES iteration.
c intype = 0 (steady)
c intype = 1 (explicit)
c intype = -1 (implicit)
include 'SIZE'
include 'TOTAL'
include 'GMRES'
common /ctolpr/ divex
common /cprint/ ifprint
logical ifprint
real res (lx2*ly2*lz2*lelv)
real h1 (lx1,ly1,lz1,lelv)
real h2 (lx1,ly1,lz1,lelv)
real h2inv(lx1,ly1,lz1,lelv)
common /scrmg/ wp (lx2,ly2,lz2,lelv)
common /ctmp0/ wk1(lgmres),wk2(lgmres)
common /cgmres1/ y(lgmres)
real alpha, l, temp
integer j,m
c
logical iflag
save iflag
data iflag /.false./
real norm_fac
save norm_fac
c
real*8 etime1,dnekclock
c
if(.not.iflag) then
iflag=.true.
call uzawa_gmres_split0(ml,mu,bm2,bm2inv,nx2*ny2*nz2*nelv)
norm_fac = 1./sqrt(volvm2)
endif
c
etime1 = dnekclock()
etime_p = 0.
divex = 0.
iter = 0
m = lgmres
c
call chktcg2(tolps,res,iconv)
if (param(21).gt.0.and.tolps.gt.abs(param(21)))
$ tolps = abs(param(21))
c if (param(21).lt.0) tolps = abs(param(21))
if (istep.eq.0) tolps = 1.e-4
tolpss = tolps
c
ntot2 = nx2*ny2*nz2*nelv
c
iconv = 0
call rzero(x,ntot2)
do while(iconv.eq.0.and.iter.lt.100)
if(iter.eq.0) then
! -1
call col3(r,ml,res,ntot2) ! r = L res
c call copy(r,res,ntot2)
else
!update residual
call copy(r,res,ntot2) ! r = res
call cdabdtp(w,x,h1,h2,h2inv,intype) ! w = A x
call add2s2(r,w,-1.,ntot2) ! r = r - w
! -1
call col2(r,ml,ntot2) ! r = L r
endif
! ______
gamma(1) = sqrt(glsc2(r,r,ntot2)) ! gamma = \/ (r,r)
! 1
if(iter.eq.0) then
div0 = gamma(1)*norm_fac
if (param(21).lt.0) tolpss=abs(param(21))*div0
endif
!check for lucky convergence
rnorm = 0.
if(gamma(1) .eq. 0.) goto 9000
temp = 1./gamma(1)
call cmult2(v(1,1),r,temp,ntot2) ! v = r / gamma
! 1 1
do j=1,m
iter = iter+1
! -1
call col3(w,mu,v(1,j),ntot2) ! w = U v
! j
etime2 = dnekclock()
if(param(43).eq.1) then
call uzprec(z(1,j),w,h1,h2,intype,wp)
else ! -1
call hsmg_solve(z(1,j),w) ! z = M w
endif
etime_p = etime_p + dnekclock()-etime2
call cdabdtp(w,z(1,j), ! w = A z
$ h1,h2,h2inv,intype) ! j
! -1
call col2(w,ml,ntot2) ! w = L w
c !modified Gram-Schmidt
c do i=1,j
c h(i,j)=glsc2(w,v(1,i),ntot2) ! h = (w,v )
c ! i,j i
c call add2s2(w,v(1,i),-h(i,j),ntot2)! w = w - h v
c enddo ! i,j i
c 2-PASS GS, 1st pass:
do i=1,j
h(i,j)=vlsc2(w,v(1,i),ntot2) ! h = (w,v )
enddo ! i,j i
call gop(h(1,j),wk1,'+ ',j) ! sum over P procs
do i=1,j
call add2s2(w,v(1,i),-h(i,j),ntot2)! w = w - h v
enddo ! i,j i
c 2-PASS GS, 2nd pass:
c
c do i=1,j
c wk1(i)=vlsc2(w,v(1,i),ntot2) ! h = (w,v )
c enddo ! i,j i
c !
c call gop(wk1,wk2,'+ ',j) ! sum over P procs
c
c do i=1,j
c call add2s2(w,v(1,i),-wk1(i),ntot2)! w = w - h v
c h(i,j) = h(i,j) + wk1(i) ! i,j i
c enddo
!apply Givens rotations to new column
do i=1,j-1
temp = h(i,j)
h(i ,j)= c(i)*temp + s(i)*h(i+1,j)
h(i+1,j)= -s(i)*temp + c(i)*h(i+1,j)
enddo
! ______
alpha = sqrt(glsc2(w,w,ntot2)) ! alpha = \/ (w,w)
rnorm = 0.
if(alpha.eq.0.) goto 900 !converged
l = sqrt(h(j,j)*h(j,j)+alpha*alpha)
temp = 1./l
c(j) = h(j,j) * temp
s(j) = alpha * temp
h(j,j) = l
gamma(j+1) = -s(j) * gamma(j)
gamma(j) = c(j) * gamma(j)
c call outmat(h,m,j,' h ',j)
rnorm = abs(gamma(j+1))*norm_fac
ratio = rnorm/div0
if (ifprint.and.nid.eq.0)
$ write (6,66) iter,tolpss,rnorm,div0,ratio,istep
66 format(i5,1p4e12.5,i8,' Divergence')
#ifndef TST_WSCAL
if (rnorm .lt. tolpss) goto 900 !converged
#else
if (iter.gt.param(151)-1) goto 900
#endif
if (j.eq.m) goto 1000 !not converged, restart
temp = 1./alpha
call cmult2(v(1,j+1),w,temp,ntot2) ! v = w / alpha
! j+1
enddo
900 iconv = 1
1000 continue
!back substitution
! -1
!c = H gamma
do k=j,1,-1
temp = gamma(k)
do i=j,k+1,-1
temp = temp - h(k,i)*c(i)
enddo
c(k) = temp/h(k,k)
enddo
!sum up Arnoldi vectors
do i=1,j
call add2s2(x,z(1,i),c(i),ntot2) ! x = x + c z
! i i
enddo
c if(iconv.eq.1) call dbg_write(x,nx2,ny2,nz2,nelv,'esol',3)
enddo
9000 continue
c
divex = rnorm
c iter = iter - 1
c
c DIAGNOSTICS
c call copy(w,x,ntot2)
c if (ifvcor) then
c xaver = glsc2(bm2,w,ntot2)/volvm2
c call cadd(w,-xaver,ntot2)
c endif
c call copy(r,res,ntot2) !r = res
c call cdabdtp(r,w,h1,h2,h2inv,intype) ! r = A w
c do i=1,ntot2
c r(i) = res(i) - r(i) ! r = res - r
c enddo
c call uzawa_gmres_temp(r,bm2inv,ntot2)
c ! ______
c gamma(1) = sqrt(glsc2(r,r,ntot2)/volvm2) ! gamma = \/ (r,r)
c ! 1
c print *, 'GMRES end resid:',gamma(1)
c END DIAGNOSTICS
call copy(res,x,ntot2)
c
if (ifvcor) then
xaver = glsc2(bm2,res,ntot2)/volvm2
call cadd(res,-xaver,ntot2)
endif
c
etime1 = dnekclock()-etime1
if (nid.eq.0) write(6,9999) istep,iter,divex,tolpss,div0,etime_p,
& etime1
c call flush_hack
9999 format(4X,I7,' U-Pres gmres: ',I6,1p5E13.4)
c
c
return
end
c-----------------------------------------------------------------------
subroutine uzawa_gmres_split0(l,u,b,binv,n)
integer n
real l(n),u(n),b(n),binv(n)
integer i
do i=1,n
l(i)=sqrt(binv(i))
u(i)=sqrt(b(i))
if(abs(u(i)*l(i)-1.0).gt.1e-13) print *, i, u(i)*l(i)
enddo
return
end
c-----------------------------------------------------------------------
subroutine uzawa_gmres_split(l,u,b,binv,n)
integer n
real l(n),u(n),b(n),binv(n)
integer i
do i=1,n
c l(i)=sqrt(binv(i))
c u(i)=sqrt(b(i))
c u(i)=sqrt(b(i))
c l(i)=1./u(i)
c l(i)=sqrt(binv(i))
l(i)=1.
u(i)=1./l(i)
c if(abs(u(i)*l(i)-1.0).gt.1e-13)write(6,*) i,u(i)*l(i),' gmr_sp'
enddo
return
end
c-----------------------------------------------------------------------
subroutine uzawa_gmres_temp(a,b,n)
integer n
real a(n),b(n)
integer i
do i=1,n
a(i)=sqrt(b(i))*a(i)
enddo
return
end
c-----------------------------------------------------------------------
subroutine ax(w,x,h1,h2,n)
include 'SIZE'
include 'TOTAL'
c
c w = A*x for pressure iteration
c
integer n
real w(n),x(n),h1(n),h2(n)
imsh = 1
isd = 1
call axhelm (w,x,h1,h2,imsh,isd)
call dssum (w,nx1,ny1,nz1)
call col2 (w,pmask,n)
return
end
c-----------------------------------------------------------------------
subroutine hmh_gmres(res,h1,h2,wt,iter)
c Solve the Helmholtz equation by right-preconditioned
c GMRES iteration.
include 'SIZE'
include 'TOTAL'
include 'FDMH1'
include 'GMRES'
common /ctolpr/ divex
common /cprint/ ifprint
logical ifprint
real res (lx1*ly1*lz1*lelv)
real h1 (lx1,ly1,lz1,lelv)
real h2 (lx1,ly1,lz1,lelv)
real wt (lx1,ly1,lz1,lelv)
common /scrcg/ d(lx1*ly1*lz1*lelv),wk(lx1*ly1*lz1*lelv)
common /cgmres1/ y(lgmres)
common /ctmp0/ wk1(lgmres),wk2(lgmres)
real alpha, l, temp
integer j,m
logical iflag
save iflag
data iflag /.false./
real norm_fac
save norm_fac
real*8 etime1,dnekclock
n = nx1*ny1*nz1*nelv
etime1 = dnekclock()
etime_p = 0.
divex = 0.
iter = 0
m = lgmres
if(.not.iflag) then
iflag=.true.
call uzawa_gmres_split(ml,mu,bm1,binvm1,nx1*ny1*nz1*nelv)
norm_fac = 1./sqrt(volvm1)
endif
call set_fdm_prec_h1b(d,h1,h2,nelv)
if (ifvcor) smean = -1./glsum(vmult,n)
call chktcg1(tolps,res,h1,h2,pmask,vmult,1,1)
if (param(21).gt.0.and.tolps.gt.abs(param(21)))
$ tolps = abs(param(21))
if (istep.eq.0) tolps = 1.e-4
tolpss = tolps
c
iconv = 0
call rzero(x,n)
do while(iconv.eq.0.and.iter.lt.500)
if(iter.eq.0) then ! -1
call col3(r,ml,res,n) ! r = L res
c call copy(r,res,n)
else
!update residual
call copy (r,res,n) ! r = res
call ax (w,x,h1,h2,n) ! w = A x
call add2s2(r,w,-1.,n) ! r = r - w
! -1
call col2(r,ml,n) ! r = L r
endif
! ______
gamma(1) = sqrt(glsc3(r,r,wt,n)) ! gamma = \/ (r,r)
! 1
if(iter.eq.0) then
div0 = gamma(1)*norm_fac
if (param(21).lt.0) tolpss=abs(param(21))*div0
endif
!check for lucky convergence
rnorm = 0.
if(gamma(1) .eq. 0.) goto 9000
temp = 1./gamma(1)
call cmult2(v(1,1),r,temp,n) ! v = r / gamma
! 1 1
do j=1,m
iter = iter+1
! -1
call col3(w,mu,v(1,j),n) ! w = U v
! j
c . . . . . Overlapping Schwarz + coarse-grid . . . . . . .
etime2 = dnekclock()
kfldfdm = ndim+1
call fdm_h1
$ (z(1,j),w,d,pmask,vmult,nelv,ktype(1,1,kfldfdm),wk)
call crs_solve_h1 (wk,w)
c call hsmg_solve(z(1,j),w) ! z = M w
! j
call add2 (z(1,j),wk,n)
if (ifvcor) then
rmean = smean*glsc2(z(1,j),vmult,n)
call cadd(z(1,j),rmean,n)
endif
etime_p = etime_p + dnekclock()-etime2
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
call ax (w,z(1,j),h1,h2,n) ! w = A z
! j
! -1
call col2(w,ml,n) ! w = L w
c !modified Gram-Schmidt
c do i=1,j
c h(i,j)=glsc3(w,v(1,i),wt,n) ! h = (w,v )
c ! i,j i
c call add2s2(w,v(1,i),-h(i,j),n) ! w = w - h v
c enddo ! i,j i
c 2-PASS GS, 1st pass:
do i=1,j
h(i,j)=vlsc3(w,v(1,i),wt,n) ! h = (w,v )
enddo ! i,j i
call gop(h(1,j),wk1,'+ ',j) ! sum over P procs
do i=1,j
call add2s2(w,v(1,i),-h(i,j),n) ! w = w - h v
enddo ! i,j i
c 2-PASS GS, 2nd pass:
c
c do i=1,j
c wk1(i)=vlsc3(w,v(1,i),wt,n) ! h = (w,v )
c enddo ! i,j i
c !
c call gop(wk1,wk2,'+ ',j) ! sum over P procs
c
c do i=1,j
c call add2s2(w,v(1,i),-wk1(i),n) ! w = w - h v
c h(i,j) = h(i,j) + wk1(i) ! i,j i
c enddo
!apply Givens rotations to new column
do i=1,j-1
temp = h(i,j)
h(i ,j)= c(i)*temp + s(i)*h(i+1,j)
h(i+1,j)= -s(i)*temp + c(i)*h(i+1,j)
enddo
! ______
alpha = sqrt(glsc3(w,w,wt,n)) ! alpha = \/ (w,w)
rnorm = 0.
if(alpha.eq.0.) goto 900 !converged
l = sqrt(h(j,j)*h(j,j)+alpha*alpha)
temp = 1./l
c(j) = h(j,j) * temp
s(j) = alpha * temp
h(j,j) = l
gamma(j+1) = -s(j) * gamma(j)
gamma(j) = c(j) * gamma(j)
rnorm = abs(gamma(j+1))*norm_fac
ratio = rnorm/div0
if (ifprint.and.nid.eq.0)
$ write (6,66) iter,tolpss,rnorm,div0,ratio,istep
66 format(i5,1p4e12.5,i8,' Divergence')
#ifndef TST_WSCAL
if (rnorm .lt. tolpss) goto 900 !converged
#else
if (iter.gt.param(151)-1) goto 900
#endif
if (j.eq.m) goto 1000 !not converged, restart
temp = 1./alpha
call cmult2(v(1,j+1),w,temp,n) ! v = w / alpha
! j+1
enddo
900 iconv = 1
1000 continue
!back substitution
! -1
!c = H gamma
do k=j,1,-1
temp = gamma(k)
do i=j,k+1,-1
temp = temp - h(k,i)*c(i)
enddo
c(k) = temp/h(k,k)
enddo
!sum up Arnoldi vectors
do i=1,j
call add2s2(x,z(1,i),c(i),n) ! x = x + c z
enddo ! i i
c if(iconv.eq.1) call dbg_write(x,nx1,ny1,nz1,nelv,'esol',3)
enddo
9000 continue
c
divex = rnorm
call copy(res,x,n)
c
if (ifvcor) then
xaver = glsc2(bm1,res,n)/volvm1
call cadd(res,-xaver,n)
endif
c
etime1 = dnekclock()-etime1
if (nid.eq.0) write(6,9999) istep,iter,divex,tolpss,div0,etime_p,
& etime1
c call flush_hack
9999 format(4X,I7,' PRES gmres: ',I6,1p5E13.4)
return
end
c-----------------------------------------------------------------------
More information about the Nek5000-users
mailing list