[Nek5000-users] Residuals for post-processing

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Sun Apr 25 12:28:49 CDT 2010


Hi,

uh-oh, I think I misread your first response - your suggestion was to 
include:
call outpost(x,y,z,res,t,'res')
call exitti('quit in gmres$',n)
in gmres.f, not in usrchk ?

What do x,v,and z stand for?

Thanks,
Markus


nek5000-users at lists.mcs.anl.gov wrote:
> 
> 
> 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
>>
> _______________________________________________
> 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