[Nek5000-users] Residuals for post-processing
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Sun Apr 25 14:21:29 CDT 2010
Hi Markus,
outpost simply provides a window towards generating
.fld files.
It's looking for 3 fields that will be interpreted as
components of the vector field (vx,vy,vz) -- but obviously
that field will contain only what is actually passed to it.
4th field is the pressure surrogate
5th field is the temperature surrogate.
It's important if using pn-pn-2 that the pressure field
look like a pressure field (the velocity and temp are pn,
pressure is pn-2).
Paul
On Sun, 25 Apr 2010, nek5000-users at lists.mcs.anl.gov wrote:
> 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
>>
> _______________________________________________
> 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