[Nek5000-users] Residuals for post-processing
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Sun Apr 25 12:49:53 CDT 2010
Hi Markus,
The issue here is that "res" is an argument passed into
the gmres solver and is not the stored value that is held
in GMRES, so a seg fault would be a natural failure.
The stored values are x,r,v,and z.
Paul
On Sun, 25 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:
> Hi,
>
> this is what I added to usrchk:
> ---
> include 'GMRES'
> if (istep.gt.0) then
> call outpost(vx,vy,vz,res,t,'res')
> call exitti('quit in gmres$',n)
> endif
> ---
> Unfortunately, this causes the code to crash whenever outpost is called:
> ---
> 120 1.00000E-06 8.11690E-04 5.58247E-01 1.45400E-03 1 Divergence
> 1 U-Pres gmres: 120 8.1169E-04 1.0000E-06 5.5825E-01
> 2.8286E+00 8.0148E+00
> 1 DNORM, DIVEX 0.35002708528548815 8.11690205515358640E-004
> 1 5.0000E-03 8.3137E+00 Fluid done
> filt amp 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0125 0.0500
> filt trn 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9875 0.9500
> ./nek: line 7: 11765 Segmentation fault (core dumped) ./nek5000
>
> real 0m27.297s
> user 0m9.348s
> sys 0m1.124s
> ---
>
> I tried with pr instead of res, that works. Is it possible that res is not an
> array but a single number?
> Where are the residuals for velocity and passive scalars stored?
>
> I added the divergence procedure to the outflow, there is no back flow. The
> problem is not that the simulation blows up (at least for the time steps I am
> doing), it is that the pressure residuals don't fall below the specified
> limit when I increase lx1. This is for both PN/PN-2 and PN/PN and not a
> function of dt.
>
> Thanks,
> Markus
>
>
>
> nek5000-users at lists.mcs.anl.gov wrote:
>>
>> 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-----------------------------------------------------------------------
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
> _______________________________________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.gov
> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>
More information about the Nek5000-users
mailing list