[Nek5000-users] Residuals for post-processing

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Wed Apr 21 10:06:42 CDT 2010


Do you end up sucking any flow into the outflow ?

If so, the std fix for that is to set div U > 0 near the
outflow...

Typically I find it sufficient to dump at every step over
the last few steps prior to blow up to identify the problem...

hth

Paul


On Wed, 21 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:

> Hi,
>
> thank you, calling it after the time step is what I need.
>
> There is a backward facing step within my domain. In order to shorten 
> reattachment length, we tried to apply suction at the backward facing wall 
> (using a velocity boundary condition with negative (opposing the mean flow 
> direction) velocity magnitude).
> Doing initial testing with PN/PN-2, this seems to work for lower polynomial 
> orders, but as soon as I up lx1 (from 6 to 8, say), gmres will not decrease 
> over iterations but remain at a constant level (this level becoming higher 
> with lx1 increasing further).
> I tried to decrease dt, but that doesn't help. So I wanted to see if high 
> pressure residuals are really in that area of the suction slot or if the 
> problem is somewhere else.
>
> Thanks again,
> 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