Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
hmholtz.F90
Go to the documentation of this file.
1 !=======================================================================
2 subroutine hmholtz(name,u,rhs,h1,h2,mask,mult,imsh,tli,maxit,isd)
3  use kinds, only : dp
4  use size_m, only : lx1, ly1, lz1, nx1, ny1, nz1, nelv, nelt, ndim, nid
5  use ctimer, only : icalld, thmhz, nhmhz, dnekclock
6  use fdmh1, only : kfldfdm
7  use input, only : ifsplit, param
8  use geom, only : binvm1, bintm1
9  use tstep, only : istep, nelfld, ifield
10  implicit none
11 
12  CHARACTER(4) NAME
13  REAL(DP) :: U (lx1,ly1,lz1,*)
14  REAL(DP) :: RHS (lx1,ly1,lz1,*)
15  REAL(DP) :: H1 (lx1,ly1,lz1,*)
16  REAL(DP) :: H2 (lx1,ly1,lz1,*)
17  REAL(DP) :: MASK (lx1,ly1,lz1,*)
18  REAL(DP) :: MULT (lx1,ly1,lz1,*)
19  real(DP) :: tli
20  integer :: imsh, maxit, isd
21 
22  logical :: iffdm
23  character(3) :: nam3
24  integer :: ntot
25  real(DP) :: tol
26  real(DP) :: etime
27 
28  tol = abs(tli)
29 
30  iffdm = .false.
31 ! iffdm = .true.
32  if (ifsplit) iffdm = .true.
33 
34 !max if (icalld == 0 .AND. iffdm) call set_fdm_prec_h1A
35 
36  nhmhz=nhmhz + 1
37  etime=dnekclock()
38  ntot = nx1*ny1*nz1*nelfld(ifield)
39  if (imsh == 1) ntot = nx1*ny1*nz1*nelv
40  if (imsh == 2) ntot = nx1*ny1*nz1*nelt
41 
42 ! Determine which field is being computed for FDM based preconditioner bc's
43 
44  call chcopy(nam3,name,3)
45 
46  kfldfdm = -1
47 ! if (nam3.eq.'TEM' ) kfldfdm = 0
48 ! if (name.eq.'TEM1') kfldfdm = 0 ! hardcode for temp only, for mpaul
49 ! if (name.eq.'VELX') kfldfdm = 1
50 ! if (name.eq.'VELY') kfldfdm = 2
51 ! if (name.eq.'VELZ') kfldfdm = 3
52  if (name == 'PRES') kfldfdm = ndim+1
53 ! if (.not.iffdm) kfldfdm=-1
54 
55  call dssum(rhs)
56  rhs(:,:,:,1:nelfld(ifield)) = rhs(:,:,:,1:nelfld(ifield)) * mask(:,:,:,1:nelfld(ifield))
57  if (nid == 0 .AND. istep <= 10) &
58  write(6,*) param(22),' p22 ',istep,imsh
59  if (param(22) == 0 .OR. istep <= 10) &
60  call chktcg1(tol,rhs,h1,h2,mask,mult,imsh,isd)
61 
62  if (tli < 0) tol=tli ! caller-specified relative tolerance
63 
64  if (imsh == 1) call cggo &
65  (u,rhs,h1,h2,mask,mult,imsh,tol,maxit,isd,binvm1,name)
66  if (imsh == 2) call cggo &
67  (u,rhs,h1,h2,mask,mult,imsh,tol,maxit,isd,bintm1,name)
68 
69 
70  thmhz=thmhz+(dnekclock()-etime)
71  return
72 end subroutine hmholtz
73 
74 !=======================================================================
75 !------------------------------------------------------------------
78 !------------------------------------------------------------------
79 subroutine axhelm (au,u,helm1,helm2,imesh,isd)
80  use kinds, only : dp
81  use size_m, only : lx1, ly1, lz1
82  use size_m, only : nx1, ny1, nz1, nelt, nelv, ndim
83  use ctimer, only : taxhm, naxhm, etime1, dnekclock
84  use ctimer, only : axhelm_flop, axhelm_mop
85  use geom, only : g4m1, g5m1, g6m1
86  use dxyz, only : wddx, wddyt, wddzt
87  use input, only : ifaxis
88  use geom, only : bm1
89  use mesh, only : ifsolv
90  implicit none
91 
92  real(DP), intent(out) :: AU (lx1,ly1,lz1,*) !>
93  real(DP), intent(in) :: U (lx1,ly1,lz1,*) !>
94  real(DP), intent(in) :: HELM1 (lx1,ly1,lz1,*) !>
95  real(DP), intent(in) :: HELM2 (lx1,ly1,lz1,*) !>
96  integer, intent(in) :: imesh !>
97  integer, intent(in) :: isd !>
98 
99  ! locals
100  REAL(DP) :: TM1 (lx1,ly1,lz1)
101  REAL(DP) :: TM2 (lx1,ly1,lz1)
102  REAL(DP) :: TM3 (lx1,ly1,lz1)
103 
104  integer :: nel, nxy, nyz, nxz, nxyz, ntot
105  integer :: e, iz
106 
107 #ifdef BGQ
108  integer :: iy, ix
109  vector(real(dp)) :: tm1v, tm2v, tm3v, g4mv, g5mv, g6mv, uv, auv, bm1v, h1v, h2v
110 #endif
111 
112  nel=nelt
113  if (imesh == 1) nel=nelv
114 
115  nxy=nx1*ny1
116  nyz=ny1*nz1
117  nxz=nx1*nz1
118  nxyz=nx1*ny1*nz1
119  ntot=nxyz*nel
120 
121  if (naxhm == 0) taxhm=0.0
122  naxhm=naxhm + 1
123 
124  IF ( .NOT. ifsolv) CALL setfast(helm1,helm2,imesh)
125  etime1=dnekclock()
126  !au(:,:,:,1:nel) = 0._dp
127 
128  do 100 e=1,nel
129 
130 ! if (ifaxis) call setaxdy ( ifrzer(e) )
131 
132  IF (ndim == 2) THEN
133  write(*,*) "Whoops! axhelm"
134 #if 0
135 
136  ! 2-d case ...............
137 
138  if (iffast(e)) then
139 
140  ! Fast 2-d mode: constant properties and undeformed element
141 
142  h1 = helm1(1,1,1,e)
143  call mxm(wddx,nx1,u(1,1,1,e),nx1,tm1,nyz)
144  call mxm(u(1,1,1,e),nx1,wddyt,ny1,tm2,ny1)
145  call col2(tm1,g4m1(1,1,1,e),nxyz)
146  call col2(tm2,g5m1(1,1,1,e),nxyz)
147  call add3(au(1,1,1,e),tm1,tm2,nxyz)
148  call cmult(au(1,1,1,e),h1,nxyz)
149 
150  else
151 
152  ! General case, speed-up for undeformed elements
153 
154  call mxm(dxm1,nx1,u(1,1,1,e),nx1,dudr,nyz)
155  call mxm(u(1,1,1,e),nx1,dytm1,ny1,duds,ny1)
156  call col3(tmp1,dudr,g1m1(1,1,1,e),nxyz)
157  call col3(tmp2,duds,g2m1(1,1,1,e),nxyz)
158  if (ifdfrm(e)) then
159  call addcol3(tmp1,duds,g4m1(1,1,1,e),nxyz)
160  call addcol3(tmp2,dudr,g4m1(1,1,1,e),nxyz)
161  endif
162  call col2(tmp1,helm1(1,1,1,e),nxyz)
163  call col2(tmp2,helm1(1,1,1,e),nxyz)
164  call mxm(dxtm1,nx1,tmp1,nx1,tm1,nyz)
165  call mxm(tmp2,nx1,dym1,ny1,tm2,ny1)
166  call add2(au(1,1,1,e),tm1,nxyz)
167  call add2(au(1,1,1,e),tm2,nxyz)
168 
169  endif
170 
171 #endif
172  else
173 
174  ! 3-d case ...............
175 
176 ! if (iffast(e)) then
177  if (.true.) then
178 
179  ! Fast 3-d mode: constant properties and undeformed element
180  axhelm_flop = axhelm_flop + (2*nx1-1)*nx1*nyz
181  axhelm_mop = axhelm_mop + nx1*ny1*nz1
182  axhelm_flop = axhelm_flop + 9*nx1*ny1*nz1
183  axhelm_mop = axhelm_mop + 5*nx1*ny1*nz1
184  axhelm_flop = axhelm_flop + (2*ny1-1)*nx1*ny1*nz1
185  axhelm_flop = axhelm_flop + nxy*(2*nz1-1)*nz1
186 
187 
188 
189  call mxm(wddx,nx1,u(1,1,1,e),nx1,tm1,nyz)
190  do iz=1,nz1
191  call mxm(u(1,1,iz,e),nx1,wddyt,ny1,tm2(1,1,iz),ny1)
192  END DO
193  call mxm(u(1,1,1,e),nxy,wddzt,nz1,tm3,nz1)
194 #ifdef BGQ
195  h1v = vec_splats(helm1(1,1,1,e))
196  h2v = vec_splats(helm2(1,1,1,e))
197  do iz = 1, nz1
198  do iy = 1, ny1
199  g4mv = vec_ld(0, g4m1(1,iy,iz,e))
200  tm1v = vec_ld(0, tm1(1,iy,iz))
201  tm1v = vec_mul(g4mv,tm1v)
202 
203  g5mv = vec_ld(0, g5m1(1,iy,iz,e))
204  tm2v = vec_ld(0, tm2(1,iy,iz))
205  tm1v = vec_madd(g5mv,tm2v,tm1v)
206 
207  g6mv = vec_ld(0, g6m1(1,iy,iz,e))
208  tm3v = vec_ld(0, tm3(1,iy,iz))
209  tm1v = vec_madd(g6mv,tm3v,tm1v)
210 
211  uv = vec_ld(0, u(1,iy,iz,e))
212  bm1v = vec_ld(0, bm1(1,iy,iz,e))
213  uv = vec_mul(uv,bm1v)
214  uv = vec_mul(uv,h2v)
215  auv = vec_madd(h1v, tm1v, uv)
216  call vec_st(auv,0,au(1,iy,iz,e))
217 
218  g4mv = vec_ld(0, g4m1(5,iy,iz,e))
219  tm1v = vec_ld(0, tm1(5,iy,iz))
220  tm1v = vec_mul(g4mv,tm1v)
221 
222  g5mv = vec_ld(0, g5m1(5,iy,iz,e))
223  tm2v = vec_ld(0, tm2(5,iy,iz))
224  tm1v = vec_madd(g5mv,tm2v,tm1v)
225 
226  g6mv = vec_ld(0, g6m1(5,iy,iz,e))
227  tm3v = vec_ld(0, tm3(5,iy,iz))
228  tm1v = vec_madd(g6mv,tm3v,tm1v)
229 
230  uv = vec_ld(0, u(5,iy,iz,e))
231  bm1v = vec_ld(0, bm1(5,iy,iz,e))
232  uv = vec_mul(uv,bm1v)
233  uv = vec_mul(uv,h2v)
234  auv = vec_madd(h1v, tm1v, uv)
235  call vec_st(auv,0,au(5,iy,iz,e))
236  enddo
237  enddo
238 #else
239  au(:,:,:,e) = helm1(1,1,1,e) * ( &
240  tm1*g4m1(:,:,:,e) &
241  + tm2*g5m1(:,:,:,e) &
242  + tm3*g6m1(:,:,:,e) ) &
243  + helm2(1,1,1,1)*bm1(:,:,:,e)*u(:,:,:,e)
244 #endif
245  else
246  write(*,*) "Woops! axhelm"
247 #if 0
248 
249  ! General case, speed-up for undeformed elements
250 
251  call mxm(dxm1,nx1,u(1,1,1,e),nx1,dudr,nyz)
252  do 10 iz=1,nz1
253  call mxm(u(1,1,iz,e),nx1,dytm1,ny1,duds(1,1,iz),ny1)
254  10 END DO
255  call mxm(u(1,1,1,e),nxy,dztm1,nz1,dudt,nz1)
256  call col3(tmp1,dudr,g1m1(1,1,1,e),nxyz)
257  call col3(tmp2,duds,g2m1(1,1,1,e),nxyz)
258  call col3(tmp3,dudt,g3m1(1,1,1,e),nxyz)
259  if (ifdfrm(e)) then
260  call addcol3(tmp1,duds,g4m1(1,1,1,e),nxyz)
261  call addcol3(tmp1,dudt,g5m1(1,1,1,e),nxyz)
262  call addcol3(tmp2,dudr,g4m1(1,1,1,e),nxyz)
263  call addcol3(tmp2,dudt,g6m1(1,1,1,e),nxyz)
264  call addcol3(tmp3,dudr,g5m1(1,1,1,e),nxyz)
265  call addcol3(tmp3,duds,g6m1(1,1,1,e),nxyz)
266  endif
267  call col2(tmp1,helm1(1,1,1,e),nxyz)
268  call col2(tmp2,helm1(1,1,1,e),nxyz)
269  call col2(tmp3,helm1(1,1,1,e),nxyz)
270  call mxm(dxtm1,nx1,tmp1,nx1,tm1,nyz)
271  do 20 iz=1,nz1
272  call mxm(tmp2(1,1,iz),nx1,dym1,ny1,tm2(1,1,iz),ny1)
273  20 END DO
274  call mxm(tmp3,nxy,dzm1,nz1,tm3,nz1)
275  call add2(au(1,1,1,e),tm1,nxyz)
276  call add2(au(1,1,1,e),tm2,nxyz)
277  call add2(au(1,1,1,e),tm3,nxyz)
278 
279 #endif
280  endif
281 
282  endif
283 
284  100 END DO
285 
286  !au(:,:,:,1:nel) = au(:,:,:,1:nel) + helm2(1,1,1,1)*bm1(:,:,:,1:nel)*u(:,:,:,1:nel)
287 ! au(:,:,:,1:nel) = au(:,:,:,1:nel) + helm2(:,:,:,1:nel)*bm1(:,:,:,1:nel)*u(:,:,:,1:nel)
288 
289 ! If axisymmetric, add a diagonal term in the radial direction (ISD=2)
290 
291  if (ifaxis .AND. (isd == 2)) then
292  write(*,*) "Whoops! axhelm 3"
293 #if 0
294  do 200 e=1,nel
295 
296  if (ifrzer(e)) then
297  call mxm(u(1,1,1,e),nx1,datm1,ny1,duax,1)
298  call mxm(ym1(1,1,1,e),nx1,datm1,ny1,ysm1,1)
299  endif
300 
301  do 190 j=1,ny1
302  do 190 i=1,nx1
303  ! if (ym1(i,j,1,e).ne.0.) then
304  if (ifrzer(e)) then
305  term1 = 0.0
306  if(j /= 1) &
307  term1 = bm1(i,j,1,e)*u(i,j,1,e)/ym1(i,j,1,e)**2
308  term2 = wxm1(i)*wam1(1)*dam1(1,j)*duax(i) &
309  *jacm1(i,1,1,e)/ysm1(i)
310  else
311  term1 = bm1(i,j,1,e)*u(i,j,1,e)/ym1(i,j,1,e)**2
312  term2 = 0.
313  endif
314  au(i,j,1,e) = au(i,j,1,e) &
315  + helm1(i,j,1,e)*(term1+term2)
316  ! endif
317  190 END DO
318  200 END DO
319 #endif
320  endif
321 
322  taxhm=taxhm+(dnekclock()-etime1)
323  return
324 end subroutine axhelm
325 
326 !=======================================================================
328 !-------------------------------------------------------------------
329 subroutine setfast (helm1,helm2,imesh)
330  use kinds, only : dp
331  use size_m, only : nx1, ny1, nz1, nelt, nelv
332  use input, only : ifaxis, ifmodel
333  use mesh, only : iffast, ifdfrm
334  use ctimer, only : nsetfast, tsetfast, dnekclock
335  implicit none
336 
337  integer, intent(in) :: imesh
338  REAL(DP), intent(in) :: HELM1(nx1,ny1,nz1,*), HELM2(nx1,ny1,nz1,*)
339 
340  integer :: nel, nxyz, ntot, ie
341  real(DP) :: delta, x, y, diff, epsm, h1min, h1max, testh1
342  real(DP) :: den
343  real(DP) :: etime
344  real(DP), external :: vlamax
345 
346  etime = dnekclock()
347  nsetfast = nsetfast + 1
348  nel = -1
349  IF (imesh == 1) nel=nelv
350  IF (imesh == 2) nel=nelt
351  nxyz = nx1*ny1*nz1
352  ntot = nxyz*nel
353 
354  delta = 1.e-9
355  x = 1.+delta
356  y = 1.
357  diff = abs(x-y)
358  IF (diff == 0.0) epsm = 1.e-6
359  IF (diff > 0.0) epsm = 1.e-13
360 
361 
362  DO 100 ie=1,nel
363  iffast(ie) = .false.
364  IF (ifdfrm(ie) .OR. ifaxis .OR. ifmodel ) THEN
365  iffast(ie) = .false.
366  ! write(*,*) IFDFRM(ie), IFAXIS, IFMODEL
367  ELSE
368  h1min = minval(helm1(:,:,:,ie))
369  h1max = maxval(helm1(:,:,:,ie))
370  den = abs(h1max)+abs(h1min)
371  if (den > 0) then
372  testh1 = abs((h1max-h1min)/(h1max+h1min))
373  IF (testh1 < epsm) iffast(ie) = .true.
374  else
375  iffast(ie) = .true.
376  endif
377  ENDIF
378  100 END DO
379  tsetfast = tsetfast + (dnekclock() - etime)
380  return
381 end subroutine setfast
382 
383 !----------------------------------------------------------------------
386 !----------------------------------------------------------------------
387 subroutine sfastax()
388  use kinds, only : dp
389  use size_m, only : nx1, ny1, nz1, nelt, ndim
390  use dxyz, only : dxm1, dym1, dzm1
391  use dxyz, only : wddx, wddyt, wddzt
392  use geom, only : g1m1, g2m1, g3m1, g4m1, g5m1, g6m1
393  use mesh, only : ifdfrm
394  use wz_m, only : wxm1, wym1, wzm1
395  implicit none
396 
397  logical, save :: IFIRST = .TRUE.
398  integer :: nxx, nyy, nzz
399  integer :: i, j, ip, ie, ix, iy, iz
400 
401  nxx=nx1*nx1
402  IF (ifirst) THEN
403  wddx = 0._dp
404  DO i=1,nx1
405  DO j=1,nx1
406  DO ip=1,nx1
407  wddx(i,j) = wddx(i,j) + wxm1(ip)*dxm1(ip,i)*dxm1(ip,j)
408  enddo
409  enddo
410  END DO
411  nyy=ny1*ny1
412  wddyt = 0._dp
413  DO i=1,ny1
414  DO j=1,ny1
415  DO ip=1,ny1
416  wddyt(j,i) = wddyt(j,i) + wym1(ip)*dym1(ip,i)*dym1(ip,j)
417  enddo
418  enddo
419  END DO
420  nzz=nz1*nz1
421  wddzt = 0._dp
422  DO i=1,nz1
423  DO j=1,nz1
424  DO ip=1,nz1
425  wddzt(j,i) = wddzt(j,i) + wzm1(ip)*dzm1(ip,i)*dzm1(ip,j)
426  enddo
427  enddo
428  END DO
429  ifirst= .false.
430  ENDIF
431 
432  IF (ndim == 3) THEN
433  DO 1001 ie=1,nelt
434  IF ( .NOT. ifdfrm(ie)) THEN
435  DO iz=1,nz1
436  DO iy=1,ny1
437  DO ix=1,nx1
438  g4m1(ix,iy,iz,ie)=g1m1(ix,iy,iz,ie)/wxm1(ix)
439  g5m1(ix,iy,iz,ie)=g2m1(ix,iy,iz,ie)/wym1(iy)
440  g6m1(ix,iy,iz,ie)=g3m1(ix,iy,iz,ie)/wzm1(iz)
441  enddo
442  enddo
443  END DO
444  ENDIF
445  1001 END DO
446  ELSE
447  DO 2001 ie=1,nelt
448  IF ( .NOT. ifdfrm(ie)) THEN
449  DO iy=1,ny1
450  DO ix=1,nx1
451  g4m1(ix,iy,1,ie)=g1m1(ix,iy,1,ie)/wxm1(ix)
452  g5m1(ix,iy,1,ie)=g2m1(ix,iy,1,ie)/wym1(iy)
453  enddo
454  END DO
455  ENDIF
456  2001 END DO
457  ENDIF
458  return
459  end subroutine sfastax
460 
461 !=======================================================================
463 !-------------------------------------------------------------------
464 subroutine setprec (dpcm1,helm1,helm2,imsh,isd)
465  use kinds, only : dp
466  use size_m, only : nx1, ny1, nz1, lx1, ly1, lz1, nelt, nelv, ndim
467  use dxyz, only : dxtm1, dytm1, dztm1, datm1, dam1
468  use geom, only : ifrzer, g1m1, g2m1, g3m1, ym1, jacm1
469  use input, only : ifaxis
470  use geom, only : bm1
471  use mesh, only : ifdfrm
472  use wz_m, only : wxm1, wam1
473  use ctimer, only : tdpc, ndpc, dnekclock
474  implicit none
475 
476  REAL(DP), intent(out) :: DPCM1 (lx1,ly1,lz1,*)
477  REAL(DP), intent(in) :: HELM1(nx1,ny1,nz1,*), HELM2(nx1,ny1,nz1,*)
478  integer, intent(in) :: imsh, isd
479 
480  REAL(DP) :: YSM1(ly1)
481 
482  integer :: nel, ntot, ie, iq, iz, iy, ix, j, i, iel
483  real(DP) :: term1, term2
484  real(DP) :: etime
485 
486  ndpc = ndpc + 1
487  etime = dnekclock()
488 
489  nel=nelt
490  if (imsh == 1) nel=nelv
491 
492  ntot = nel*nx1*ny1*nz1
493 
494 ! The following lines provide a convenient debugging option
495 ! call rone(dpcm1,ntot)
496 ! if (ifield.eq.1) call copy(dpcm1,binvm1,ntot)
497 ! if (ifield.eq.2) call copy(dpcm1,bintm1,ntot)
498 ! return
499 
500  dpcm1(:,:,:,1:nel) = 0._dp
501  DO 1000 ie=1,nel
502 
503 ! IF (IFAXIS) CALL SETAXDY ( IFRZER(IE) )
504 
505  DO iz=1,nz1
506  DO iy=1,ny1
507  DO ix=1,nx1
508  DO iq=1,nx1
509  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) + &
510  g1m1(iq,iy,iz,ie) * dxtm1(ix,iq)**2
511  END DO
512  enddo
513  enddo
514  enddo
515  DO iz=1,nz1
516  DO iy=1,ny1
517  DO ix=1,nx1
518  DO iq=1,ny1
519  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) + &
520  g2m1(ix,iq,iz,ie) * dytm1(iy,iq)**2
521  END DO
522  enddo
523  enddo
524  enddo
525 
526  IF (ndim == 3) THEN
527  DO iz=1,nz1
528  DO iy=1,ny1
529  DO ix=1,nx1
530  DO iq=1,nz1
531  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) + &
532  g3m1(ix,iy,iq,ie) * dztm1(iz,iq)**2
533  END DO
534  enddo
535  enddo
536  enddo
537 
538  ! Add cross terms if element is deformed.
539 
540  IF (ifdfrm(ie)) THEN
541  write(*,*) "Whoops!"
542 #if 0
543  DO 600 iy=1,ny1
544  DO 600 iz=1,nz1
545  dpcm1(1,iy,iz,ie) = dpcm1(1,iy,iz,ie) &
546  + g4m1(1,iy,iz,ie) * dxtm1(1,1)*dytm1(iy,iy) &
547  + g5m1(1,iy,iz,ie) * dxtm1(1,1)*dztm1(iz,iz)
548  dpcm1(nx1,iy,iz,ie) = dpcm1(nx1,iy,iz,ie) &
549  + g4m1(nx1,iy,iz,ie) * dxtm1(nx1,nx1)*dytm1(iy,iy) &
550  + g5m1(nx1,iy,iz,ie) * dxtm1(nx1,nx1)*dztm1(iz,iz)
551  600 END DO
552  DO 700 ix=1,nx1
553  DO 700 iz=1,nz1
554  dpcm1(ix,1,iz,ie) = dpcm1(ix,1,iz,ie) &
555  + g4m1(ix,1,iz,ie) * dytm1(1,1)*dxtm1(ix,ix) &
556  + g6m1(ix,1,iz,ie) * dytm1(1,1)*dztm1(iz,iz)
557  dpcm1(ix,ny1,iz,ie) = dpcm1(ix,ny1,iz,ie) &
558  + g4m1(ix,ny1,iz,ie) * dytm1(ny1,ny1)*dxtm1(ix,ix) &
559  + g6m1(ix,ny1,iz,ie) * dytm1(ny1,ny1)*dztm1(iz,iz)
560  700 END DO
561  DO 800 ix=1,nx1
562  DO 800 iy=1,ny1
563  dpcm1(ix,iy,1,ie) = dpcm1(ix,iy,1,ie) &
564  + g5m1(ix,iy,1,ie) * dztm1(1,1)*dxtm1(ix,ix) &
565  + g6m1(ix,iy,1,ie) * dztm1(1,1)*dytm1(iy,iy)
566  dpcm1(ix,iy,nz1,ie) = dpcm1(ix,iy,nz1,ie) &
567  + g5m1(ix,iy,nz1,ie) * dztm1(nz1,nz1)*dxtm1(ix,ix) &
568  + g6m1(ix,iy,nz1,ie) * dztm1(nz1,nz1)*dytm1(iy,iy)
569  800 END DO
570 #endif
571  ENDIF
572  ELSE
573 
574  IF (ifdfrm(ie)) THEN
575  write(*,*) "Whoops!"
576 #if 0
577  iz=1
578  DO 602 iy=1,ny1
579  dpcm1(1,iy,iz,ie) = dpcm1(1,iy,iz,ie) &
580  + g4m1(1,iy,iz,ie) * dxtm1(1,1)*dytm1(iy,iy)
581  dpcm1(nx1,iy,iz,ie) = dpcm1(nx1,iy,iz,ie) &
582  + g4m1(nx1,iy,iz,ie) * dxtm1(nx1,nx1)*dytm1(iy,iy)
583  602 END DO
584  DO 702 ix=1,nx1
585  DO 702 iz=1,nz1
586  dpcm1(ix,1,iz,ie) = dpcm1(ix,1,iz,ie) &
587  + g4m1(ix,1,iz,ie) * dytm1(1,1)*dxtm1(ix,ix)
588  dpcm1(ix,ny1,iz,ie) = dpcm1(ix,ny1,iz,ie) &
589  + g4m1(ix,ny1,iz,ie) * dytm1(ny1,ny1)*dxtm1(ix,ix)
590  702 END DO
591 #endif
592  ENDIF
593  ENDIF
594  1000 END DO
595 
596  dpcm1(:,:,:,1:nel) = dpcm1(:,:,:,1:nel)*helm1(:,:,:,1:nel) + bm1*helm2(:,:,:,1:nel)
597 
598 ! If axisymmetric, add a diagonal term in the radial direction (ISD=2)
599 
600  IF (ifaxis .AND. (isd == 2)) THEN
601  DO 1200 iel=1,nel
602 
603  IF (ifrzer(iel)) THEN
604  CALL mxm(ym1(1,1,1,iel),nx1,datm1,ny1,ysm1,1)
605  ENDIF
606 
607  DO j=1,ny1
608  DO i=1,nx1
609  IF (ym1(i,j,1,iel) /= 0.) THEN
610  term1 = bm1(i,j,1,iel)/ym1(i,j,1,iel)**2
611  IF (ifrzer(iel)) THEN
612  term2 = wxm1(i)*wam1(1)*dam1(1,j) &
613  *jacm1(i,1,1,iel)/ysm1(i)
614  ELSE
615  term2 = 0.
616  ENDIF
617  dpcm1(i,j,1,iel) = dpcm1(i,j,1,iel) &
618  + helm1(i,j,1,iel)*(term1+term2)
619  ENDIF
620  enddo
621  END DO
622  1200 END DO
623  ENDIF
624 
625  CALL dssum(dpcm1)
626  dpcm1(:,:,:,1:nel) = 1._dp / dpcm1(:,:,:,1:nel)
627  tdpc = tdpc + (dnekclock() - etime)
628 
629  return
630 end subroutine setprec
631 
632 !-------------------------------------------------------------------
636 !-------------------------------------------------------------------
637 subroutine chktcg1 (tol,res,h1,h2,mask,mult,imesh,isd)
638  use kinds, only : dp
639  use size_m, only : lx1, ly1, lz1
640  use size_m, only : nx1, ny1, nz1, nelv, nelt, nid
641  use eigen, only : eigaa, eigga
642  use input, only : ifprint
643  use geom, only : volvm1, voltm1, binvm1, bintm1, bm1
644  use ctimer, only : nfoo, tfoo, dnekclock
645  implicit none
646 
647  real(DP), intent(out) :: tol
648  REAL(DP), intent(in) :: RES (lx1,ly1,lz1,*)
649  REAL(DP), intent(in) :: H1 (lx1,ly1,lz1,*)
650  REAL(DP), intent(in) :: H2 (lx1,ly1,lz1,*)
651  REAL(DP), intent(in) :: MASK (lx1,ly1,lz1,*)
652  REAL(DP), intent(in) :: MULT (lx1,ly1,lz1,*)
653  integer, intent(in) :: imesh, isd
654 
655  real(DP), allocatable :: W1(:,:,:,:), W2(:,:,:,:)
656 
657  real(DP) :: acondno, delta, x, y, diff, eps
658  real(DP) :: vol, rinit, rmin, bcneu1, bcneu2, bctest, bcrob, tolmin
659  real(DP) :: etime
660  integer :: nl, ntot1
661  real(DP), external :: glsc3, glsum
662 
663  nfoo = nfoo + 1
664  etime = dnekclock()
665 
666  IF (eigaa /= 0.) THEN
667  acondno = eigga/eigaa
668  ELSE
669  acondno = 10.
670  ENDIF
671 
672 ! Single or double precision???
673 
674  delta = 1.e-9
675  x = 1.+delta
676  y = 1.
677  diff = abs(x-y)
678  IF (diff == 0.) eps = 1.e-6
679  IF (diff > 0.) eps = 1.e-13
680  ! max
681  IF (diff > 0.) eps = 1.e-16
682 
683  IF (imesh == 1) THEN
684  nl = nelv
685  vol = volvm1
686  ELSEIF (imesh == 2) THEN
687  nl = nelt
688  vol = voltm1
689  else
690  nl = -1
691  vol = -1.0
692  write(*,*) "Got somewhere bad"
693  ENDIF
694 
695  allocate(w1(nx1,ny1,nz1,nl), w2(nx1,ny1,nz1,nl))
696  ntot1 = nx1*ny1*nz1*nl
697 
698  IF (imesh == 1) THEN
699  w2 = binvm1 * res(:,:,:,1:nl)
700  rinit = sqrt(glsc3(w2,res,mult,ntot1)/volvm1)
701  ELSE
702  w2 = bintm1 * res(:,:,:,1:nl)
703  rinit = sqrt(glsc3(w2,res,mult,ntot1)/voltm1)
704  ENDIF
705  rmin = eps*rinit
706  IF (tol < rmin) THEN
707  IF (nid == 0 .AND. ifprint) &
708  WRITE (6,*) 'New CG1-tolerance (RINIT*epsm) = ',rmin,tol
709  tol = rmin
710  ENDIF
711 
712  w1 = 1._dp
713  bcneu1 = glsc3(w1,mask,mult,ntot1)
714  bcneu2 = glsc3(w1,w1 ,mult,ntot1)
715  bctest = abs(bcneu1-bcneu2)
716 
717  etime = etime - dnekclock()
718  CALL axhelm(w2,w1,h1,h2,imesh,isd)
719  etime = etime + dnekclock()
720  w2 = w2 * w2 * bm1
721  bcrob = sqrt(glsum(w2,ntot1)/vol)
722 
723  IF ((bctest < .1) .AND. (bcrob < (eps*acondno))) THEN
724  ! OTR = GLSC3 (W1,RES,MULT,NTOT1)
725  tolmin = rinit*eps*10.
726  IF (tol < tolmin) THEN
727  tol = tolmin
728  IF (nid == 0 .AND. ifprint) &
729  WRITE(6,*) 'New CG1-tolerance (Neumann) = ',tolmin
730  ENDIF
731  ENDIF
732 
733  tfoo = tfoo + (dnekclock() - etime)
734 
735  return
736  end subroutine chktcg1
737 
738 !=======================================================================
739 !-------------------------------------------------------------------------
743 !------------------------------------------------------------------------
744 subroutine cggo(x,f,h1,h2,mask,mult,imsh,tin,maxit,isd,binv,name)
745  use kinds, only : dp
746  use size_m, only : nid, nx1, ny1, nz1, nelt, nelv
747  use size_m, only : lx1, ly1, lz1, lelt
748  use fdmh1, only : kfldfdm
749  use input, only : ifsplit, param, ifprint
750  use geom, only : volvm1, voltm1
751  use mesh, only : ifsolv, niterhm
752  use tstep, only : istep, imesh
753  use ctimer, only : ncggo, tcggo, cggo_flop, cggo_mop, dnekclock
754  implicit none
755 
756  real(DP), intent(out) :: x(lx1*ly1*lz1,*) !>
757  real(DP), intent(in) :: f(lx1*ly1*lz1,*) !>
758  real(DP), intent(in) :: h1(lx1*ly1*lz1,*) !>
759  real(DP), intent(in) :: h2(lx1*ly1*lz1,*) !>
760  real(DP), intent(in) :: mask(lx1*ly1*lz1,*) !>
761  real(DP), intent(in) :: mult(lx1*ly1*lz1,*) !>
762  real(DP), intent(in) :: binv(lx1*ly1*lz1,*) !>
763  integer, intent(in) :: imsh, isd, maxit
764  real(DP), intent(in) :: tin !>
765  character(4) :: name
766 
767  logical :: ifmcor,ifprint_hmh
768 
769  integer, parameter :: lxyz = lx1*ly1*lz1
770  real(DP) :: scalar(2)
771  real(DP), allocatable :: d (:,:)
772  real(DP), allocatable :: r (:,:) , w (:,:) , p (:,:) , z (:,:)
773 
774  integer :: i, n, iter, nxyz, nel, niter
775  real(DP) :: rho, vol, tol, h2max, skmin
776  real(DP) :: rtz1, rtz2, rbn2, rbn0, beta, rho0, alpha, alphm
777  real(DP), external :: glmax, glmin, glsum, glsc2, vlsc3, vlsc32, glsc3
778  real(DP) :: etime
779 
780  if (ifsplit .AND. name == 'PRES' .AND. param(42) == 0) then
781  n = nx1*ny1*nz1*nelv
782  call copy(x,f,n)
783  call hmh_gmres(x,h1,h2,mult,iter)
784  niterhm = iter
785  return
786  endif
787 ! write(6,*) ifsplit,name,param(44),' P44 C'
788 
789 ! ** zero out stuff for Lanczos eigenvalue estimator
790  ncggo = ncggo + 1
791 
792  rho = 0.00
793 
794 ! Initialization
795  nxyz = nx1*ny1*nz1
796  nel = nelv
797  vol = volvm1
798  IF (imsh == 2) nel=nelt
799  IF (imsh == 2) vol=voltm1
800  n = nel*nxyz
801 
802  tol=abs(tin)
803  if (param(22) /= 0) tol=abs(param(22))
804  if (name == 'PRES' .AND. param(21) /= 0) tol=abs(param(21))
805  if (tin < 0) tol=abs(tin)
806  niter = maxit
807 
808 ! Speed-up for undeformed elements and constant properties.
809  if ( .NOT. ifsolv) then
810  call setfast(h1,h2,imesh)
811  ifsolv = .true.
812  endif
813 
814 ! Set up diag preconditioner.
815 
816  allocate(d(lxyz,lelt)); d = 0_dp
817  if (kfldfdm < 0) then
818  call setprec(d,h1,h2,imsh,isd)
819  elseif(param(100) /= 2) then
820 !max call set_fdm_prec_h1b(d,h1,h2,nel)
821  endif
822  etime = dnekclock()
823 
824  allocate(r(nxyz,nel))
825  cggo_mop = cggo_mop + 3*n
826  r = f(:,1:nel)
827  x(:,1:nel) = 0._dp
828 
829 ! Check for non-trivial null-space
830 
831  ifmcor = .false.
832  cggo_mop = cggo_mop + 2*n
833  h2max = glmax(h2 ,n)
834  skmin = glmin(mask,n)
835  if (skmin > 0 .AND. h2max == 0) then
836  ifmcor = .true.
837  write(*,*) "Oops! ifmcor"
838  endif
839  if (name == 'PRES') then
840  write(*,*) "Oops! name == PRES"
841  endif
842  if (kfldfdm >= 0) then
843  write(*,*) "Oops! kfldfdm"
844  endif
845 
846  ifprint_hmh = .false.
847  if (nid == 0 .AND. ifprint .AND. param(74) /= 0) ifprint_hmh= .true.
848  if (nid == 0 .AND. istep == 1) ifprint_hmh= .true.
849 
850  niterhm = 0
851 
852  cggo_mop = cggo_mop + n
853  allocate(w(nxyz,nel)); w = 0_dp
854  allocate(p(nxyz,nel))
855  iter = 1
856 
857  rtz1=1._dp; rtz2 = 1._dp
858 
859  cggo_flop = cggo_flop + 8*n
860  cggo_mop = cggo_mop + 5*n
861  scalar(1) = 0._dp; scalar(2) = 0._dp
862  do i = 1, nel
863  p(:,i) = r(:,i) * d(:,i)
864  scalar(1) = scalar(1) + sum(p(:,i) * r(:,i) * mult(:,i))
865  scalar(2) = scalar(2) + sum(r(:,i) * r(:,i) * mult(:,i) * binv(:,i))
866  enddo
867  call gop(scalar,w,'+ ',2)
868  rtz1=scalar(1)
869  rbn2=sqrt(scalar(2)/vol)
870  rbn0 = rbn2
871  if (param(22) < 0) tol=abs(param(22))*rbn0
872  if (tin < 0) tol=abs(tin)*rbn0
873 
874  if (ifprint_hmh) &
875  write(6,3002) istep,iter,name,ifmcor,rbn2,tol,h1(1,1),h2(1,1)
876 
877  beta = 0._dp
878 
879  etime = etime - dnekclock()
880  call axhelm(w,p,h1,h2,imsh,isd)
881  etime = etime + dnekclock()
882  call dssum(w)
883 
884  cggo_flop = cggo_flop + 4*n
885  cggo_mop = cggo_mop + 5*n
886  rho0 = rho; rho = 0._dp
887  do i = 1, nel
888  w(:,i) = w(:,i) * mask(:,i)
889  rho = rho + sum(w(:,i) * p(:,i) * mult(:,i))
890  enddo
891  allocate(z(nxyz,nel))
892  call gop(rho,z,'+ ',1)
893  alpha=rtz1/rho
894  alphm=-alpha
895 
896  do iter=2,niter
897  scalar(1) = 0._dp; scalar(2) = 0._dp
898  rtz2=rtz1
899  cggo_flop = cggo_flop + 10*n
900  cggo_mop = cggo_mop + 7*n
901  do i = 1, nel
902  r(:,i) = r(:,i) + alphm * w(:,i)
903  z(:,i) = r(:,i) * d(:,i)
904  scalar(1) = scalar(1) + sum(z(:,i) * r(:,i) * mult(:,i))
905  scalar(2) = scalar(2) + sum(r(:,i) * r(:,i) * mult(:,i) * binv(:,i))
906  enddo
907  call gop(scalar,w,'+ ',2)
908  rtz1=scalar(1)
909  rbn2=sqrt(scalar(2)/vol)
910 
911  if (ifprint_hmh) &
912  write(6,3002) istep,iter,name,ifmcor,rbn2,tol,h1(1,1),h2(1,1)
913 
914 
915  ! Always take at least one iteration (for projection) pff 11/23/98
916 #ifndef TST_WSCAL
917  !IF (rbn2 <= TOL .AND. iter > 10) THEN
918  IF (rbn2 <= tol .AND. (iter > 1 .OR. istep <= 5)) THEN
919 #else
920  iter_max = param(150)
921  if (name == 'PRES') iter_max = param(151)
922  if (iter > iter_max) then
923 #endif
924  ! IF (rbn2.LE.TOL) THEN
925  cggo_flop = cggo_flop + 2*n
926  cggo_mop = cggo_mop + 3*n
927  x(:,1:nel) = x(:,1:nel) + alpha * p
928  niter = iter-1
929  ! IF(NID.EQ.0.AND.((.NOT.IFHZPC).OR.IFPRINT))
930  if (nid == 0) &
931  write(6,3000) istep,name,niter,rbn2,rbn0,tol
932  goto 9999
933  ENDIF
934 
935  cggo_flop = cggo_flop + 4*n
936  cggo_mop = cggo_mop + 5*n
937  beta = rtz1/rtz2
938  do i = 1, nel
939  x(:,i) = x(:,i) + alpha * p(:,i)
940  p(:,i) = z(:,i) + beta * p(:,i)
941  enddo
942  etime = etime - dnekclock()
943  call axhelm(w,p,h1,h2,imsh,isd)
944  etime = etime + dnekclock()
945  call dssum(w)
946 
947  cggo_flop = cggo_flop + 4*n
948  cggo_mop = cggo_mop + 5*n
949  rho0 = rho; rho = 0._dp
950  do i = 1, nel
951  w(:,i) = w(:,i) * mask(:,i)
952  rho = rho + sum(w(:,i) * p(:,i) * mult(:,i))
953  enddo
954  call gop(rho,z,'+ ',1)
955  alpha=rtz1/rho
956  alphm=-alpha
957  enddo
958 
959  niter = iter-1
960 
961  if (nid == 0) write (6,3001) istep,niter,name,rbn2,rbn0,tol
962  3000 format(4x,i7,4x,'Hmholtz ',a4,': ',i6,1p6e13.4)
963  3001 format(2i6,' **ERROR**: Failed in HMHOLTZ: ',a4,1p6e13.4)
964  3002 format(i3,i6,' Helmholtz ',a4,1x,l4,':',1p6e13.4)
965  9999 continue
966  niterhm = niter
967  ifsolv = .false.
968  tcggo = tcggo + (dnekclock() - etime)
969 
970  return
971 end subroutine cggo
972 
973 !=======================================================================
974 real(DP) function vlsc32(r,b,m,n)
975  use kinds, only : dp
976  implicit none
977  integer, intent(in) :: n
978  real(DP), intent(in) :: r(n),b(n),m(n)
979 
980  integer :: i
981  vlsc32 = 0.
982  do i=1,n
983  vlsc32 = vlsc32 + b(i)*m(i)*r(i)*r(i)
984  enddo
985  return
986 end function vlsc32
987 
988 !=======================================================================
989 subroutine set_fdm_prec_h1b(d,h1,h2,nel)
990  use kinds, only : dp
991  use size_m, only : nx1, ny1, nz1
992  implicit none
993 
994  real(DP), intent(out) :: d (nx1,ny1,nz1,1)
995  real(DP), intent(in) :: h1(nx1,ny1,nz1,1)
996  real(DP), intent(in) :: h2(nx1,ny1,nz1,1)
997  integer, intent(in) :: nel
998 
999  d = 0._dp
1000  return
1001 
1002 #if 0
1003 ! use kinds, only : DP, DP_eps
1004 ! use size_m, only : nx1, ny1, nz1
1005 ! use fdmh1, only : ktype, elsize, dd, ifbhalf, kfldfdm
1006 ! use geom, only : ym1
1007 ! use input, only : if3d, ifaxis
1008 
1009 
1010  integer :: nxyz, i1, i2, ie, i3, k1, k2, k3
1011  real(DP) :: h1b, h2b, vol, vl1, vl2, vl3, den
1012  real(DP), external :: vlsum, vlsc2
1013 
1014 
1015 ! Set up diagonal for FDM for each spectral element
1016  nxyz = nx1*ny1*nz1
1017  if (if3d) then
1018  do ie=1,nel
1019  h1b = vlsum(h1(1,1,1,ie),nxyz)/nxyz
1020  h2b = vlsum(h2(1,1,1,ie),nxyz)/nxyz
1021  k1 = ktype(ie,1,kfldfdm)
1022  k2 = ktype(ie,2,kfldfdm)
1023  k3 = ktype(ie,3,kfldfdm)
1024  vol = elsize(1,ie)*elsize(2,ie)*elsize(3,ie)
1025  vl1 = elsize(2,ie)*elsize(3,ie)/(elsize(1,ie)+dp_eps)
1026  vl2 = elsize(1,ie)*elsize(3,ie)/(elsize(2,ie)+dp_eps)
1027  vl3 = elsize(1,ie)*elsize(2,ie)/(elsize(3,ie)+dp_eps)
1028  do i3=1,nz1
1029  do i2=1,ny1
1030  do i1=1,nx1
1031  den = h1b*(vl1*dd(i1,k1) + vl2*dd(i2,k2) + vl3*dd(i3,k3)) &
1032  + h2b*vol
1033  if (ifbhalf) den = den/vol
1034  if (den /= 0) then
1035  d(i1,i2,i3,ie) = 1./den
1036  else
1037  d(i1,i2,i3,ie) = 0.
1038 
1039  ! write(6,3) 'd=0:'
1040  ! $ ,h1(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2),dd(i3,k3)
1041  ! $ ,i1,i2,i3,ie,kfldfdm,k1,k2,k3
1042  3 format(a4,1p4e12.4,8i8)
1043 
1044  endif
1045  enddo
1046  enddo
1047  enddo
1048  enddo
1049  else
1050  do ie=1,nel
1051  if (ifaxis) then
1052  h1b = vlsc2(h1(1,1,1,ie),ym1(1,1,1,ie),nxyz)/nxyz
1053  h2b = vlsc2(h2(1,1,1,ie),ym1(1,1,1,ie),nxyz)/nxyz
1054  else
1055  h1b = vlsum(h1(1,1,1,ie),nxyz)/nxyz
1056  h2b = vlsum(h2(1,1,1,ie),nxyz)/nxyz
1057  endif
1058  k1 = ktype(ie,1,kfldfdm)
1059  k2 = ktype(ie,2,kfldfdm)
1060  vol = elsize(1,ie)*elsize(2,ie)
1061  vl1 = elsize(2,ie)/elsize(1,ie)
1062  vl2 = elsize(1,ie)/elsize(2,ie)
1063  i3=1
1064  do i2=1,ny1
1065  do i1=1,nx1
1066  den = h1b*( vl1*dd(i1,k1) + vl2*dd(i2,k2) ) &
1067  + h2b*vol
1068  if (ifbhalf) den = den/vol
1069  if (den /= 0) then
1070  d(i1,i2,i3,ie) = 1./den
1071  ! write(6,3) 'dn0:'
1072  ! $ ,d(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2)
1073  ! $ ,i1,i2,i3,ie,kfldfdm,k1,k2
1074  else
1075  d(i1,i2,i3,ie) = 0.
1076  ! write(6,3) 'd=0:'
1077  ! $ ,h1(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2)
1078  ! $ ,i1,i2,i3,ie,kfldfdm,k1,k2
1079  2 format(a4,1p3e12.4,8i8)
1080  endif
1081  ! write(6,1) ie,i1,i2,k1,k2,'d:',d(i1,i2,i3,ie),vol,vl1,vl2
1082  ! 1 format(5i3,2x,a2,1p4e12.4)
1083  enddo
1084  enddo
1085  enddo
1086  endif
1087 #endif
1088 
1089  return
1090 end subroutine set_fdm_prec_h1b
1091 
1092 !-----------------------------------------------------------------------
1101 subroutine generalev(a,b,lam,n,w)
1102  use kinds, only : dp
1103  use size_m, only : lx1, ly1, lz1, lelv, nid
1104  use parallel, only : ifdblas
1105  implicit none
1106 
1107  integer, intent(in) :: n
1108  real(DP), intent(inout) :: a(n,n), b(n,n)
1109  real(DP), intent(out) :: lam(n)
1110  real(DP), intent(in) :: w(n,n) !>
1111 
1112  real(DP) :: aa(100),bb(100)
1113 
1114  integer, parameter :: lbw=4*lx1*ly1*lz1*lelv
1115  real(DP), allocatable :: bw(:)
1116 
1117  integer :: info, ninf, lw
1118 
1119  allocate(bw(lbw))
1120 
1121  lw = n*n
1122 ! write(6,*) 'in generalev, =',info,n,ninf
1123 
1124  call copy(aa,a,100)
1125  call copy(bb,b,100)
1126 
1127  if (ifdblas) then
1128  call dsygv(1,'V','U',n,a,n,b,n,lam,bw,lbw,info)
1129  else
1130  call ssygv(1,'V','U',n,a,n,b,n,lam,bw,lbw,info)
1131  endif
1132 
1133  if (info /= 0) then
1134  if (nid == 0) then
1135  call outmat2(aa ,n,n,n,'aa ')
1136  call outmat2(bb ,n,n,n,'bb ')
1137  call outmat2(a ,n,n,n,'Aeig')
1138  call outmat2(lam,1,n,n,'Deig')
1139  endif
1140 
1141  ninf = n-info
1142  write(6,*) 'Error in generalev, info=',info,n,ninf
1143  call exitt
1144  endif
1145 
1146  return
1147 end subroutine generalev
1148 
1149 !-----------------------------------------------------------------------
1150 subroutine outmat2(a,m,n,k,name)
1151  use size_m, only : nid
1152  use kinds, only : dp
1153  implicit none
1154  integer :: m,n,k
1155  real(DP) :: a(m,n)
1156  character(4) :: name
1157  integer :: n2, i, j
1158 
1159  n2 = min(n,8)
1160  write(6,2) nid,name,m,n,k
1161  do i=1,m
1162  write(6,1) nid,name,(a(i,j),j=1,n2)
1163  enddo
1164 ! 1 format(i3,1x,a4,16f6.2)
1165  1 format(i3,1x,a4,1p8e14.5)
1166  2 format(/,'Matrix: ',i3,1x,a4,3i8)
1167  return
1168 end subroutine outmat2
1169 
1170 !-----------------------------------------------------------------------
sum
static double sum(struct xxt *data, double v, uint n, uint tag)
Definition: xxt.c:400
set_fdm_prec_h1b
subroutine set_fdm_prec_h1b(d, h1, h2, nel)
Definition: hmholtz.F90:989
outmat2
subroutine outmat2(a, m, n, k, name)
Definition: hmholtz.F90:1150
dssum
subroutine dssum(u)
Direct stiffness sum.
Definition: dssum.F90:54
vlsc32
real(dp) function vlsc32(r, b, m, n)
Definition: hmholtz.F90:974
tstep
cleaned
Definition: tstep_mod.F90:2
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
wz_m
cleaned
Definition: wz_mod.F90:2
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product C = A*B.
Definition: mxm_wrapper.F90:7
dxyz
Definition: dxyz_mod.F90:3
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Compute the (Helmholtz) matrix-vector product, AU = helm1*[A]u + helm2*[B]u, for NEL elements...
Definition: hmholtz.F90:79
exitt
void exitt()
Definition: comm_mpi.F90:411
mesh
cleaned
Definition: mesh_mod.F90:2
ctimer::dnekclock
real(dp) function dnekclock()
Definition: ctimer_mod.F90:103
copy
subroutine copy(a, b, n)
Definition: math.F90:52
sfastax
subroutine sfastax()
For undeformed elements, set up appropriate elemental matrices and geometric factors for fast evaluat...
Definition: hmholtz.F90:387
setprec
subroutine setprec(dpcm1, helm1, helm2, imsh, isd)
Generate diagonal preconditioner for the Helmholtz operator.
Definition: hmholtz.F90:464
generalev
subroutine generalev(a, b, lam, n, w)
Solve the generalized eigenvalue problem A x = lam B x.
Definition: hmholtz.F90:1101
hmh_gmres
subroutine hmh_gmres(res, h1, h2, wt, iter)
Solve the Helmholtz equation by right-preconditioned GMRES iteration.
Definition: gmres.F90:27
chktcg1
subroutine chktcg1(tol, res, h1, h2, mask, mult, imesh, isd)
Check that the tolerances are not too small for the CG-solver. Important when calling the CG-solver (...
Definition: hmholtz.F90:637
parallel
cleaned
Definition: parallel_mod.F90:2
i
for i
Definition: xxt_test.m:74
setfast
subroutine setfast(helm1, helm2, imesh)
Set logicals for fast evaluation of A*x.
Definition: hmholtz.F90:329
geom
Geometry arrays.
Definition: geom_mod.F90:2
gop
subroutine gop(x, w, op, n)
Global vector commutative operation.
Definition: comm_mpi.F90:104
cggo
subroutine cggo(x, f, h1, h2, mask, mult, imsh, tin, maxit, isd, binv, name)
Solve the Helmholtz equation, H*U = RHS, using preconditioned conjugate gradient iteration. Preconditioner: diag(H).
Definition: hmholtz.F90:744
chcopy
subroutine chcopy(a, b, n)
Definition: math.F90:63
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.F90:2
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8