[Nek5000-users] Residuals for post-processing
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Sun Apr 25 10:05:09 CDT 2010
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
>
More information about the Nek5000-users
mailing list