KTH framework for Nek5000 toolboxes; testing version  0.0.1
hmholtz.f
Go to the documentation of this file.
1 c=======================================================================
2  subroutine hmholtz(name,u,rhs,h1,h2,mask,mult,imsh,tli,maxit,isd)
3  include 'SIZE'
4  include 'TOTAL'
5  include 'FDMH1'
6  include 'CTIMER'
7 
8  CHARACTER NAME*4
9  REAL U (LX1,LY1,LZ1,1)
10  REAL RHS (LX1,LY1,LZ1,1)
11  REAL H1 (LX1,LY1,LZ1,1)
12  REAL H2 (LX1,LY1,LZ1,1)
13  REAL MASK (LX1,LY1,LZ1,1)
14  REAL MULT (LX1,LY1,LZ1,1)
15 
16  logical iffdm
17  character*3 nam3
18 
19  tol = abs(tli)
20 
21  iffdm = .false.
22  if (ifsplit) iffdm = .true.
23  if (icalld.eq.0.and.iffdm) call set_fdm_prec_h1a
24  icalld = icalld+1
25 
26 #ifdef TIMER
27  if (name.ne.'PRES') then
28  nhmhz = nhmhz + 1
29  etime1 = dnekclock()
30  endif
31 #endif
32 
33  ntot = lx1*ly1*lz1*nelfld(ifield)
34  if (imsh.eq.1) ntot = lx1*ly1*lz1*nelv
35  if (imsh.eq.2) ntot = lx1*ly1*lz1*nelt
36 
37 C Determine which field is being computed for FDM based preconditioner bc's
38 c
39  call chcopy(nam3,name,3)
40 c
41  kfldfdm = -1
42 c if (nam3.eq.'TEM' ) kfldfdm = 0
43 c if (name.eq.'TEM1') kfldfdm = 0 ! hardcode for temp only, for mpaul
44 c if (name.eq.'VELX') kfldfdm = 1
45 c if (name.eq.'VELY') kfldfdm = 2
46 c if (name.eq.'VELZ') kfldfdm = 3
47  if (name.eq.'PRES') kfldfdm = ldim+1
48 c if (.not.iffdm) kfldfdm=-1
49 C
50  call dssum (rhs,lx1,ly1,lz1)
51  call col2 (rhs,mask,ntot)
52 c if (nio.eq.0.and.istep.le.10)
53 c $ write(6,*) param(22),' p22 ',istep,imsh
54  if (param(22).eq.0.or.istep.le.10)
55  $ call chktcg1 (tol,rhs,h1,h2,mask,mult,imsh,isd)
56 
57  if (tli.lt.0) tol=tli ! caller-specified relative tolerance
58 
59  if (imsh.eq.1) call cggo
60  $ (u,rhs,h1,h2,mask,mult,imsh,tol,maxit,isd,binvm1,name)
61  if (imsh.eq.2) call cggo
62  $ (u,rhs,h1,h2,mask,mult,imsh,tol,maxit,isd,bintm1,name)
63 
64 #ifdef TIMER
65  if (name.ne.'PRES') thmhz=thmhz+(dnekclock()-etime1)
66 #endif
67 
68  return
69  END
70 C
71 c=======================================================================
72  subroutine axhelm (au,u,helm1,helm2,imesh,isd)
73 C------------------------------------------------------------------
74 C
75 C Compute the (Helmholtz) matrix-vector product,
76 C AU = helm1*[A]u + helm2*[B]u, for NEL elements.
77 C
78 C------------------------------------------------------------------
79  include 'SIZE'
80  include 'WZ'
81  include 'DXYZ'
82  include 'GEOM'
83  include 'MASS'
84  include 'INPUT'
85  include 'PARALLEL'
86  include 'CTIMER'
87 C
88  COMMON /fastax/ wddx(lx1,lx1),wddyt(ly1,ly1),wddzt(lz1,lz1)
89  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
90  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
91 C
92  REAL AU (LX1,LY1,LZ1,1)
93  $ , U (LX1,LY1,LZ1,1)
94  $ , HELM1 (LX1,LY1,LZ1,1)
95  $ , HELM2 (LX1,LY1,LZ1,1)
96  COMMON /ctmp1/ dudr(lx1,ly1,lz1)
97  $ , duds(lx1,ly1,lz1)
98  $ , dudt(lx1,ly1,lz1)
99  $ , tmp1(lx1,ly1,lz1)
100  $ , tmp2(lx1,ly1,lz1)
101  $ , tmp3(lx1,ly1,lz1)
102 
103  REAL TM1 (LX1,LY1,LZ1)
104  REAL TM2 (LX1,LY1,LZ1)
105  REAL TM3 (LX1,LY1,LZ1)
106  REAL DUAX (LX1)
107  REAL YSM1 (LX1)
108  equivalence(dudr,tm1),(duds,tm2),(dudt,tm3)
109 
110  integer e
111 
112  naxhm = naxhm + 1
113  etime1 = dnekclock()
114 
115  nel=nelt
116  if (imesh.eq.1) nel=nelv
117 
118  nxy=lx1*ly1
119  nyz=ly1*lz1
120  nxz=lx1*lz1
121  nxyz=lx1*ly1*lz1
122  ntot=nxyz*nel
123 
124  IF (.NOT.ifsolv) CALL setfast(helm1,helm2,imesh)
125  CALL rzero (au,ntot)
126 
127  do 100 e=1,nel
128 C
129  if (ifaxis) call setaxdy ( ifrzer(e) )
130 C
131  IF (ldim.EQ.2) THEN
132 C
133 C 2-d case ...............
134 C
135  if (iffast(e)) then
136 C
137 C Fast 2-d mode: constant properties and undeformed element
138 C
139  h1 = helm1(1,1,1,e)
140  call mxm (wddx,lx1,u(1,1,1,e),lx1,tm1,nyz)
141  call mxm (u(1,1,1,e),lx1,wddyt,ly1,tm2,ly1)
142  call col2 (tm1,g4m1(1,1,1,e),nxyz)
143  call col2 (tm2,g5m1(1,1,1,e),nxyz)
144  call add3 (au(1,1,1,e),tm1,tm2,nxyz)
145  call cmult (au(1,1,1,e),h1,nxyz)
146 C
147  else
148 C
149 C
150  call mxm (dxm1,lx1,u(1,1,1,e),lx1,dudr,nyz)
151  call mxm (u(1,1,1,e),lx1,dytm1,ly1,duds,ly1)
152  call col3 (tmp1,dudr,g1m1(1,1,1,e),nxyz)
153  call col3 (tmp2,duds,g2m1(1,1,1,e),nxyz)
154  if (ifdfrm(e)) then
155  call addcol3 (tmp1,duds,g4m1(1,1,1,e),nxyz)
156  call addcol3 (tmp2,dudr,g4m1(1,1,1,e),nxyz)
157  endif
158  call col2 (tmp1,helm1(1,1,1,e),nxyz)
159  call col2 (tmp2,helm1(1,1,1,e),nxyz)
160  call mxm (dxtm1,lx1,tmp1,lx1,tm1,nyz)
161  call mxm (tmp2,lx1,dym1,ly1,tm2,ly1)
162  call add2 (au(1,1,1,e),tm1,nxyz)
163  call add2 (au(1,1,1,e),tm2,nxyz)
164 
165  endif
166 C
167  else
168 C
169 C 3-d case ...............
170 C
171  if (iffast(e)) then
172 C
173 C Fast 3-d mode: constant properties and undeformed element
174 C
175  h1 = helm1(1,1,1,e)
176  call mxm (wddx,lx1,u(1,1,1,e),lx1,tm1,nyz)
177  do 5 iz=1,lz1
178  call mxm (u(1,1,iz,e),lx1,wddyt,ly1,tm2(1,1,iz),ly1)
179  5 continue
180  call mxm (u(1,1,1,e),nxy,wddzt,lz1,tm3,lz1)
181  call col2 (tm1,g4m1(1,1,1,e),nxyz)
182  call col2 (tm2,g5m1(1,1,1,e),nxyz)
183  call col2 (tm3,g6m1(1,1,1,e),nxyz)
184  call add3 (au(1,1,1,e),tm1,tm2,nxyz)
185  call add2 (au(1,1,1,e),tm3,nxyz)
186  call cmult (au(1,1,1,e),h1,nxyz)
187 C
188  else
189 C
190 C
191  call mxm(dxm1,lx1,u(1,1,1,e),lx1,dudr,nyz)
192  do 10 iz=1,lz1
193  call mxm(u(1,1,iz,e),lx1,dytm1,ly1,duds(1,1,iz),ly1)
194  10 continue
195  call mxm (u(1,1,1,e),nxy,dztm1,lz1,dudt,lz1)
196  call col3 (tmp1,dudr,g1m1(1,1,1,e),nxyz)
197  call col3 (tmp2,duds,g2m1(1,1,1,e),nxyz)
198  call col3 (tmp3,dudt,g3m1(1,1,1,e),nxyz)
199  if (ifdfrm(e)) then
200  call addcol3 (tmp1,duds,g4m1(1,1,1,e),nxyz)
201  call addcol3 (tmp1,dudt,g5m1(1,1,1,e),nxyz)
202  call addcol3 (tmp2,dudr,g4m1(1,1,1,e),nxyz)
203  call addcol3 (tmp2,dudt,g6m1(1,1,1,e),nxyz)
204  call addcol3 (tmp3,dudr,g5m1(1,1,1,e),nxyz)
205  call addcol3 (tmp3,duds,g6m1(1,1,1,e),nxyz)
206  endif
207  call col2 (tmp1,helm1(1,1,1,e),nxyz)
208  call col2 (tmp2,helm1(1,1,1,e),nxyz)
209  call col2 (tmp3,helm1(1,1,1,e),nxyz)
210  call mxm (dxtm1,lx1,tmp1,lx1,tm1,nyz)
211  do 20 iz=1,lz1
212  call mxm(tmp2(1,1,iz),lx1,dym1,ly1,tm2(1,1,iz),ly1)
213  20 continue
214  call mxm (tmp3,nxy,dzm1,lz1,tm3,lz1)
215  call add2 (au(1,1,1,e),tm1,nxyz)
216  call add2 (au(1,1,1,e),tm2,nxyz)
217  call add2 (au(1,1,1,e),tm3,nxyz)
218 C
219  endif
220 c
221  endif
222 C
223  100 continue
224 C
225  if (ifh2) call addcol4 (au,helm2,bm1,u,ntot)
226 C
227 C If axisymmetric, add a diagonal term in the radial direction (ISD=2)
228 C
229  if (ifaxis.and.(isd.eq.2)) then
230  do 200 e=1,nel
231 C
232  if (ifrzer(e)) then
233  call mxm(u(1,1,1,e),lx1,datm1,ly1,duax,1)
234  call mxm(ym1(1,1,1,e),lx1,datm1,ly1,ysm1,1)
235  endif
236 c
237  do 190 j=1,ly1
238  do 190 i=1,lx1
239 C if (ym1(i,j,1,e).ne.0.) then
240  if (ifrzer(e)) then
241  term1 = 0.0
242  if(j.ne.1)
243  $ term1 = bm1(i,j,1,e)*u(i,j,1,e)/ym1(i,j,1,e)**2
244  term2 = wxm1(i)*wam1(1)*dam1(1,j)*duax(i)
245  $ *jacm1(i,1,1,e)/ysm1(i)
246  else
247  term1 = bm1(i,j,1,e)*u(i,j,1,e)/ym1(i,j,1,e)**2
248  term2 = 0.
249  endif
250  au(i,j,1,e) = au(i,j,1,e)
251  $ + helm1(i,j,1,e)*(term1+term2)
252 C endif
253  190 continue
254  200 continue
255  endif
256 
257  taxhm=taxhm+(dnekclock()-etime1)
258  return
259  end
260 C
261 c=======================================================================
262  subroutine setfast (helm1,helm2,imesh)
263 C-------------------------------------------------------------------
264 C
265 C Set logicals for fast evaluation of A*x
266 C
267 C-------------------------------------------------------------------
268  include 'SIZE'
269  include 'INPUT'
270  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
271  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
272  REAL HELM1(lx1,ly1,lz1,1), HELM2(lx1,ly1,lz1,1)
273 C
274  IF (imesh.EQ.1) nel=nelv
275  IF (imesh.EQ.2) nel=nelt
276  nxyz = lx1*ly1*lz1
277  ntot = nxyz*nel
278 C
279  delta = 1.e-9
280  x = 1.+delta
281  y = 1.
282  diff = abs(x-y)
283  IF (diff.EQ.0.0) epsm = 1.e-6
284  IF (diff.GT.0.0) epsm = 1.e-13
285 C
286  DO 100 ie=1,nel
287  iffast(ie) = .false.
288  IF (ifdfrm(ie).OR.ifaxis) THEN
289  iffast(ie) = .false.
290  ELSE
291  h1min = vlmin(helm1(1,1,1,ie),nxyz)
292  h1max = vlmax(helm1(1,1,1,ie),nxyz)
293  den = abs(h1max)+abs(h1min)
294  if (den.gt.0) then
295  testh1 = abs((h1max-h1min)/(h1max+h1min))
296  IF (testh1.LT.epsm) iffast(ie) = .true.
297  else
298  iffast(ie) = .true.
299  endif
300  ENDIF
301  100 CONTINUE
302 c
303  ifh2 = .false.
304  testh2 = vlamax(helm2,ntot)
305  IF (testh2.GT.0.) ifh2 = .true.
306  return
307  END
308 C
309 c=======================================================================
310  subroutine sfastax
311 C----------------------------------------------------------------------
312 C
313 C For undeformed elements, set up appropriate elemental matrices
314 C and geometric factors for fast evaluation of Ax.
315 C
316 C----------------------------------------------------------------------
317  include 'SIZE'
318  include 'WZ'
319  include 'DXYZ'
320  include 'GEOM'
321  COMMON /fastax/ wddx(lx1,ly1),wddyt(ly1,ly1),wddzt(lz1,lz1)
322  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
323  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
324  LOGICAL IFIRST
325  SAVE ifirst
326  DATA ifirst /.true./
327 C
328  nxx=lx1*lx1
329  IF (ifirst) THEN
330  CALL rzero(wddx,nxx)
331  DO 100 i=1,lx1
332  DO 100 j=1,lx1
333  DO 100 ip=1,lx1
334  wddx(i,j) = wddx(i,j) + wxm1(ip)*dxm1(ip,i)*dxm1(ip,j)
335  100 CONTINUE
336  nyy=ly1*ly1
337  CALL rzero(wddyt,nyy)
338  DO 200 i=1,ly1
339  DO 200 j=1,ly1
340  DO 200 ip=1,ly1
341  wddyt(j,i) = wddyt(j,i) + wym1(ip)*dym1(ip,i)*dym1(ip,j)
342  200 CONTINUE
343  nzz=lz1*lz1
344  CALL rzero(wddzt,nzz)
345  DO 300 i=1,lz1
346  DO 300 j=1,lz1
347  DO 300 ip=1,lz1
348  wddzt(j,i) = wddzt(j,i) + wzm1(ip)*dzm1(ip,i)*dzm1(ip,j)
349  300 CONTINUE
350  ifirst=.false.
351  ENDIF
352 C
353  IF (ldim.EQ.3) THEN
354  DO 1001 ie=1,nelt
355  IF (.NOT.ifdfrm(ie)) THEN
356  DO 1000 iz=1,lz1
357  DO 1000 iy=1,ly1
358  DO 1000 ix=1,lx1
359  g4m1(ix,iy,iz,ie)=g1m1(ix,iy,iz,ie)/wxm1(ix)
360  g5m1(ix,iy,iz,ie)=g2m1(ix,iy,iz,ie)/wym1(iy)
361  g6m1(ix,iy,iz,ie)=g3m1(ix,iy,iz,ie)/wzm1(iz)
362  1000 CONTINUE
363  ENDIF
364  1001 CONTINUE
365  ELSE
366  DO 2001 ie=1,nelt
367  IF (.NOT.ifdfrm(ie)) THEN
368  DO 2000 iy=1,ly1
369  DO 2000 ix=1,lx1
370  g4m1(ix,iy,1,ie)=g1m1(ix,iy,1,ie)/wxm1(ix)
371  g5m1(ix,iy,1,ie)=g2m1(ix,iy,1,ie)/wym1(iy)
372  2000 CONTINUE
373  ENDIF
374  2001 CONTINUE
375  ENDIF
376  return
377  END
378 C
379 c=======================================================================
380  subroutine setprec (dpcm1,helm1,helm2,imsh,isd)
381 C-------------------------------------------------------------------
382 C
383 C Generate diagonal preconditioner for the Helmholtz operator.
384 C
385 C-------------------------------------------------------------------
386  include 'SIZE'
387  include 'WZ'
388  include 'DXYZ'
389  include 'GEOM'
390  include 'INPUT'
391  include 'TSTEP'
392  include 'MASS'
393  REAL DPCM1 (LX1,LY1,LZ1,1)
394  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
395  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
396  REAL HELM1(lx1,ly1,lz1,1), HELM2(lx1,ly1,lz1,1)
397  REAL YSM1(LY1)
398 
399  nel=nelt
400  if (imsh.eq.1) nel=nelv
401 
402  ntot = nel*lx1*ly1*lz1
403 
404 c The following lines provide a convenient debugging option
405 c call rone(dpcm1,ntot)
406 c if (ifield.eq.1) call copy(dpcm1,binvm1,ntot)
407 c if (ifield.eq.2) call copy(dpcm1,bintm1,ntot)
408 c return
409 
410  CALL rzero(dpcm1,ntot)
411  DO 1000 ie=1,nel
412 
413  IF (ifaxis) CALL setaxdy ( ifrzer(ie) )
414 
415  DO 320 iq=1,lx1
416  DO 320 iz=1,lz1
417  DO 320 iy=1,ly1
418  DO 320 ix=1,lx1
419  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) +
420  $ g1m1(iq,iy,iz,ie) * dxtm1(ix,iq)**2
421  320 CONTINUE
422  DO 340 iq=1,ly1
423  DO 340 iz=1,lz1
424  DO 340 iy=1,ly1
425  DO 340 ix=1,lx1
426  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) +
427  $ g2m1(ix,iq,iz,ie) * dytm1(iy,iq)**2
428  340 CONTINUE
429  IF (ldim.EQ.3) THEN
430  DO 360 iq=1,lz1
431  DO 360 iz=1,lz1
432  DO 360 iy=1,ly1
433  DO 360 ix=1,lx1
434  dpcm1(ix,iy,iz,ie) = dpcm1(ix,iy,iz,ie) +
435  $ g3m1(ix,iy,iq,ie) * dztm1(iz,iq)**2
436  360 CONTINUE
437 C
438 C Add cross terms if element is deformed.
439 C
440  IF (ifdfrm(ie)) THEN
441  DO 600 iy=1,ly1,ly1-1
442  DO 600 iz=1,lz1,max(1,lz1-1)
443  dpcm1(1,iy,iz,ie) = dpcm1(1,iy,iz,ie)
444  $ + g4m1(1,iy,iz,ie) * dxtm1(1,1)*dytm1(iy,iy)
445  $ + g5m1(1,iy,iz,ie) * dxtm1(1,1)*dztm1(iz,iz)
446  dpcm1(lx1,iy,iz,ie) = dpcm1(lx1,iy,iz,ie)
447  $ + g4m1(lx1,iy,iz,ie) * dxtm1(lx1,lx1)*dytm1(iy,iy)
448  $ + g5m1(lx1,iy,iz,ie) * dxtm1(lx1,lx1)*dztm1(iz,iz)
449  600 CONTINUE
450  DO 700 ix=1,lx1,lx1-1
451  DO 700 iz=1,lz1,max(1,lz1-1)
452  dpcm1(ix,1,iz,ie) = dpcm1(ix,1,iz,ie)
453  $ + g4m1(ix,1,iz,ie) * dytm1(1,1)*dxtm1(ix,ix)
454  $ + g6m1(ix,1,iz,ie) * dytm1(1,1)*dztm1(iz,iz)
455  dpcm1(ix,ly1,iz,ie) = dpcm1(ix,ly1,iz,ie)
456  $ + g4m1(ix,ly1,iz,ie) * dytm1(ly1,ly1)*dxtm1(ix,ix)
457  $ + g6m1(ix,ly1,iz,ie) * dytm1(ly1,ly1)*dztm1(iz,iz)
458  700 CONTINUE
459  DO 800 ix=1,lx1,lx1-1
460  DO 800 iy=1,ly1,ly1-1
461  dpcm1(ix,iy,1,ie) = dpcm1(ix,iy,1,ie)
462  $ + g5m1(ix,iy,1,ie) * dztm1(1,1)*dxtm1(ix,ix)
463  $ + g6m1(ix,iy,1,ie) * dztm1(1,1)*dytm1(iy,iy)
464  dpcm1(ix,iy,lz1,ie) = dpcm1(ix,iy,lz1,ie)
465  $ + g5m1(ix,iy,lz1,ie) * dztm1(lz1,lz1)*dxtm1(ix,ix)
466  $ + g6m1(ix,iy,lz1,ie) * dztm1(lz1,lz1)*dytm1(iy,iy)
467  800 CONTINUE
468  ENDIF
469 
470  ELSE ! 2D
471 
472  iz=1
473  IF (ifdfrm(ie)) THEN
474  DO 602 iy=1,ly1,ly1-1
475  dpcm1(1,iy,iz,ie) = dpcm1(1,iy,iz,ie)
476  $ + g4m1(1,iy,iz,ie) * dxtm1(1,1)*dytm1(iy,iy)
477  dpcm1(lx1,iy,iz,ie) = dpcm1(lx1,iy,iz,ie)
478  $ + g4m1(lx1,iy,iz,ie) * dxtm1(lx1,lx1)*dytm1(iy,iy)
479  602 CONTINUE
480  DO 702 ix=1,lx1,lx1-1
481  dpcm1(ix,1,iz,ie) = dpcm1(ix,1,iz,ie)
482  $ + g4m1(ix,1,iz,ie) * dytm1(1,1)*dxtm1(ix,ix)
483  dpcm1(ix,ly1,iz,ie) = dpcm1(ix,ly1,iz,ie)
484  $ + g4m1(ix,ly1,iz,ie) * dytm1(ly1,ly1)*dxtm1(ix,ix)
485  702 CONTINUE
486  ENDIF
487 
488  ENDIF
489  1000 CONTINUE
490 C
491  CALL col2 (dpcm1,helm1,ntot)
492  CALL addcol3 (dpcm1,helm2,bm1,ntot)
493 C
494 C If axisymmetric, add a diagonal term in the radial direction (ISD=2)
495 C
496  IF (ifaxis.AND.(isd.EQ.2)) THEN
497  DO 1200 iel=1,nel
498 C
499  IF (ifrzer(iel)) THEN
500  CALL mxm(ym1(1,1,1,iel),lx1,datm1,ly1,ysm1,1)
501  ENDIF
502 C
503  DO 1190 j=1,ly1
504  DO 1190 i=1,lx1
505  IF (ym1(i,j,1,iel).NE.0.) THEN
506  term1 = bm1(i,j,1,iel)/ym1(i,j,1,iel)**2
507  IF (ifrzer(iel)) THEN
508  term2 = wxm1(i)*wam1(1)*dam1(1,j)
509  $ *jacm1(i,1,1,iel)/ysm1(i)
510  ELSE
511  term2 = 0.
512  ENDIF
513  dpcm1(i,j,1,iel) = dpcm1(i,j,1,iel)
514  $ + helm1(i,j,1,iel)*(term1+term2)
515  ENDIF
516  1190 CONTINUE
517  1200 CONTINUE
518  ENDIF
519 C
520  CALL dssum (dpcm1,lx1,ly1,lz1)
521  CALL invcol1 (dpcm1,ntot)
522 C
523  return
524  END
525 C
526 c=======================================================================
527  subroutine chktcg1 (tol,res,h1,h2,mask,mult,imesh,isd)
528 C-------------------------------------------------------------------
529 C
530 C Check that the tolerances are not too small for the CG-solver.
531 C Important when calling the CG-solver (Gauss-Lobatto mesh) with
532 C zero Neumann b.c.
533 C
534 C-------------------------------------------------------------------
535  include 'SIZE'
536  include 'INPUT'
537  include 'MASS'
538  include 'EIGEN'
539  COMMON /cprint/ ifprint
540  LOGICAL IFPRINT
541  COMMON /ctmp0/ w1(lx1,ly1,lz1,lelt)
542  $ , w2(lx1,ly1,lz1,lelt)
543  REAL RES (LX1,LY1,LZ1,1)
544  REAL H1 (LX1,LY1,LZ1,1)
545  REAL H2 (LX1,LY1,LZ1,1)
546  REAL MULT (LX1,LY1,LZ1,1)
547  REAL MASK (LX1,LY1,LZ1,1)
548 C
549  IF (eigaa.NE.0.) THEN
550  acondno = eigga/eigaa
551  ELSE
552  acondno = 10.
553  ENDIF
554 C
555 C Single or double precision???
556 C
557  delta = 1.e-9
558  x = 1.+delta
559  y = 1.
560  diff = abs(x-y)
561  IF (diff.EQ.0.) eps = 1.e-6
562  IF (diff.GT.0.) eps = 1.e-13
563 C
564  IF (imesh.EQ.1) THEN
565  nl = nelv
566  vol = volvm1
567  ELSEIF (imesh.EQ.2) THEN
568  nl = nelt
569  vol = voltm1
570  ENDIF
571  ntot1 = lx1*ly1*lz1*nl
572  CALL copy (w1,res,ntot1)
573 C
574  IF (imesh.EQ.1) THEN
575  CALL col3 (w2,binvm1,w1,ntot1)
576  rinit = sqrt(glsc3(w2,w1,mult,ntot1)/volvm1)
577  ELSE
578  CALL col3 (w2,bintm1,w1,ntot1)
579  rinit = sqrt(glsc3(w2,w1,mult,ntot1)/voltm1)
580  ENDIF
581  rmin = eps*rinit
582  IF (tol.LT.rmin) THEN
583  IF (nio.EQ.0.AND.ifprint)
584  $ WRITE (6,*) 'New CG1-tolerance (RINIT*epsm) = ',rmin,tol
585  tol = rmin
586  ENDIF
587 C
588  CALL rone (w1,ntot1)
589  bcneu1 = glsc3(w1,mask,mult,ntot1)
590  bcneu2 = glsc3(w1,w1 ,mult,ntot1)
591  bctest = abs(bcneu1-bcneu2)
592 C
593  CALL axhelm (w2,w1,h1,h2,imesh,isd)
594  CALL col2 (w2,w2,ntot1)
595  CALL col2 (w2,bm1,ntot1)
596  bcrob = sqrt(glsum(w2,ntot1)/vol)
597 C
598  IF ((bctest .LT. .1).AND.(bcrob.LT.(eps*acondno))) THEN
599 C OTR = GLSC3 (W1,RES,MULT,NTOT1)
600  tolmin = rinit*eps*10.
601  IF (tol .LT. tolmin) THEN
602  tol = tolmin
603  IF (nio.EQ.0.AND.ifprint)
604  $ WRITE(6,*) 'New CG1-tolerance (Neumann) = ',tolmin
605  ENDIF
606  ENDIF
607 C
608  return
609  end
610 c=======================================================================
611  subroutine cggo(x,f,h1,h2,mask,mult,imsh,tin,maxit,isd,binv,name)
612 C-------------------------------------------------------------------------
613 C
614 C Solve the Helmholtz equation, H*U = RHS,
615 C using preconditioned conjugate gradient iteration.
616 C Preconditioner: diag(H).
617 C
618 C------------------------------------------------------------------------
619  include 'SIZE'
620  include 'TOTAL'
621  include 'FDMH1'
622 
623  COMMON /cprint/ ifprint, ifhzpc
624  LOGICAL IFPRINT, IFHZPC
625 
626  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
627  logical ifdfrm, iffast, ifh2, ifsolv
628 
629  logical ifmcor,ifprint_hmh
630 
631  real x(1),f(1),h1(1),h2(1),mask(1),mult(1),binv(1)
632  parameter(lg=lx1*ly1*lz1*lelt)
633  COMMON /scrcg/ d(lg) , scalar(2)
634  common /scrmg/ r(lg) , w(lg) , p(lg) , z(lg)
635 c
636  parameter(maxcg=900)
637  common /tdarray/ diagt(maxcg),upper(maxcg)
638  common /iterhm/ niterhm
639  character*4 name
640 c
641  if (ifsplit.and.name.eq.'PRES') then
642  if (param(42).eq.0) then
643  n = lx1*ly1*lz1*nelv
644  call copy (x,f,n)
645  iter = maxit
646  call hmh_gmres (x,h1,h2,mult,iter)
647  niterhm = iter
648  return
649  elseif(param(42).eq.2) then
650  n = lx1*ly1*lz1*nelv
651  call copy (x,f,n)
652  iter = maxit
653  call hmh_flex_cg(x,h1,h2,mult,iter)
654  niterhm = iter
655  return
656  endif
657  endif
658 
659 c ** zero out stuff for Lanczos eigenvalue estimator
660  call rzero(diagt,maxcg)
661  call rzero(upper,maxcg)
662  rho = 0.00
663 C
664 C Initialization
665 C
666  nxyz = lx1*ly1*lz1
667  nel = nelv
668  vol = volvm1
669  IF (imsh.EQ.2) nel=nelt
670  IF (imsh.EQ.2) vol=voltm1
671  n = nel*nxyz
672 
673  tol=abs(tin)
674 
675 c overrule input tolerance
676  if (restol(ifield).ne.0) tol=restol(ifield)
677  if (name.eq.'PRES'.and.param(21).ne.0) tol=abs(param(21))
678 
679  if (tin.lt.0) tol=abs(tin)
680  niter = min(maxit,maxcg)
681 
682  if (.not.ifsolv) then
683  call setfast(h1,h2,imsh)
684  ifsolv = .true.
685  endif
686 C
687 C Set up diag preconditioner.
688 C
689  if (kfldfdm.lt.0) then
690  call setprec(d,h1,h2,imsh,isd)
691  elseif(param(100).ne.2) then
692  call set_fdm_prec_h1b(d,h1,h2,nel)
693  endif
694 
695  call copy (r,f,n)
696  call rzero(x,n)
697  call rzero(p,n)
698 
699  fmax = glamax(f,n)
700  if (fmax.eq.0.0) return
701 
702 c Check for non-trivial null-space
703 
704  ifmcor = .false.
705  h2max = glmax(h2 ,n)
706  skmin = glmin(mask,n)
707  if (skmin.gt.0.and.h2max.eq.0) ifmcor = .true.
708 C
709  if (name.eq.'PRES') then
710 c call ortho (r) ! Commented out March 15, 2011,pff
711  elseif (ifmcor) then
712 
713  smean = -1./glsum(bm1,n) ! Modified 5/4/12 pff
714  rmean = smean*glsc2(r,mult,n)
715  call copy(x,bm1,n)
716  call dssum(x,lx1,ly1,lz1)
717  call add2s2(r,x,rmean,n)
718  call rzero(x,n)
719  endif
720 C
721  krylov = 0
722  rtz1=1.0
723  niterhm = 0
724 
725  do iter=1,niter
726 C
727  if (kfldfdm.lt.0) then ! Jacobi Preconditioner
728 c call copy(z,r,n)
729  call col3(z,r,d,n)
730  else ! Schwarz Preconditioner
731  if (name.eq.'PRES'.and.param(100).eq.2) then
732  call h1_overlap_2(z,r,mask)
733  call crs_solve_h1 (w,r) ! Currently, crs grd only for P
734  call add2 (z,w,n)
735  else
736  call fdm_h1(z,r,d,mask,mult,nel,ktype(1,1,kfldfdm),w)
737  if (name.eq.'PRES') then
738  call crs_solve_h1 (w,r) ! Currently, crs grd only for P
739  call add2 (z,w,n)
740  endif
741  endif
742  endif
743 c
744  if (name.eq.'PRES') then
745  call ortho (z)
746  elseif (ifmcor) then
747  rmean = smean*glsc2(z,bm1,n)
748  call cadd(z,rmean,n)
749  endif
750 c write(6,*) rmean,ifmcor,' ifmcor'
751 c
752  rtz2=rtz1
753  scalar(1)=vlsc3(z,r,mult,n)
754  if(param(18).eq.1) then
755  scalar(2)=vlsc3(r,r,mult,n)
756  else
757  scalar(2)=vlsc32(r,mult,binv,n)
758  endif
759  call gop(scalar,w,'+ ',2)
760  rtz1=scalar(1)
761  rbn2=sqrt(scalar(2)/vol)
762  if (iter.eq.1) rbn0 = rbn2
763  if (param(22).lt.0) tol=abs(param(22))*rbn0
764  if (tin.lt.0) tol=abs(tin)*rbn0
765 
766  ifprint_hmh = .false.
767  if (nio.eq.0.and.ifprint.and.param(74).ne.0) ifprint_hmh=.true.
768  if (nio.eq.0.and.istep.eq.1) ifprint_hmh=.true.
769 
770  if (ifprint_hmh)
771  & write(6,3002) istep,' Hmholtz ' // name,
772  & iter,rbn2,h1(1),tol,h2(1),ifmcor
773 
774 
775 c Always take at least one iteration (for projection) pff 11/23/98
776 #ifndef FIXITER
777  IF (rbn2.LE.tol.and.(iter.gt.1 .or. istep.le.5)) THEN
778 #else
779  iter_max = param(150)
780  if (name.eq.'PRES') iter_max = param(151)
781  if (iter.gt.iter_max) then
782 #endif
783 c IF (rbn2.LE.TOL) THEN
784  niter = iter-1
785 c IF(NID.EQ.0.AND.((.NOT.IFHZPC).OR.IFPRINT))
786  if (nio.eq.0)
787  & write(6,3000) istep,' Hmholtz ' // name,
788  & niter,rbn2,rbn0,tol
789  goto 9999
790  ENDIF
791 c
792  beta = rtz1/rtz2
793  if (iter.eq.1) beta=0.0
794  call add2s1 (p,z,beta,n)
795  call axhelm (w,p,h1,h2,imsh,isd)
796  call dssum (w,lx1,ly1,lz1)
797  call col2 (w,mask,n)
798 c
799  rho0 = rho
800  rho = glsc3(w,p,mult,n)
801  alpha=rtz1/rho
802  alphm=-alpha
803  call add2s2(x,p ,alpha,n)
804  call add2s2(r,w ,alphm,n)
805 c
806 c Generate tridiagonal matrix for Lanczos scheme
807  if (iter.eq.1) then
808  krylov = krylov+1
809  diagt(iter) = rho/rtz1
810  elseif (iter.le.maxcg) then
811  krylov = krylov+1
812  diagt(iter) = (beta**2 * rho0 + rho ) / rtz1
813  upper(iter-1) = -beta * rho0 / sqrt(rtz2 * rtz1)
814  endif
815  1000 enddo
816  niter = iter-1
817 c
818  if (nio.eq.0) write (6,3001) istep, ' Error Hmholtz ' // name,
819  & niter,rbn2,rbn0,tol
820 
821 
822  3000 format(i11,a,1x,i7,1p4e13.4)
823  3001 format(i11,a,1x,i7,1p4e13.4)
824  3002 format(i11,a,1x,i7,1p4e13.4,l4)
825  9999 continue
826  niterhm = niter
827  ifsolv = .false.
828 c
829 c
830 c Call eigenvalue routine for Lanczos scheme:
831 c two work arrays are req'd if you want to save "diag & upper"
832 c
833 c if (iter.ge.3) then
834 c niter = iter-1
835 c call calc (diagt,upper,w,z,krylov,dmax,dmin)
836 c cond = dmax/dmin
837 c if (nid.eq.0) write(6,6) istep,cond,dmin,dmax,' lambda'
838 c endif
839 c 6 format(i9,1p3e12.4,4x,a7)
840 c
841 c if (n.gt.0) write(6,*) 'quit in cggo'
842 c if (n.gt.0) call exitt
843 c call exitt
844  return
845  end
846 c=======================================================================
847  function vlsc32(r,b,m,n)
848  real r(1),b(1),m(1)
849  s = 0.
850  do i=1,n
851  s = s + b(i)*m(i)*r(i)*r(i)
852  enddo
853  vlsc32 = s
854  return
855  end
856 c=======================================================================
857  subroutine calc (diag,upper,d,e,n,dmax,dmin)
858 c
859  dimension diag(n),upper(n)
860  dimension d(n),e(n)
861 c
862  call copy (d,diag ,n)
863  call copy (e,upper,n)
864 c
865  do 15 l=1,n
866  iter = 0
867 c
868  1 do 12 m=l,n-1
869  dd = abs( d(m) ) + abs( d(m+1) )
870  if ( abs(e(m)) + dd .eq. dd ) goto 2
871  12 continue
872 c
873  m = n
874  2 if ( m .ne. l ) then
875 c
876  if ( iter .eq. 30 ) then
877  write (6,*) 'too many iterations'
878  return
879  endif
880 c
881  iter = iter + 1
882  g = ( d(l+1) - d(l) ) / ( 2.0 * e(l) )
883  r = sqrt( g**2 + 1.0 )
884 c
885 c sign is defined as a(2) * abs( a(1) )
886 c
887  g = d(m) - d(l) + e(l)/(g+sign(r,g))
888  s = 1.0
889  c = 1.0
890  p = 0.0
891 c
892  do 14 i = m-1,l,-1
893  f = s * e(i)
894  b = c * e(i)
895  if ( abs(f) .ge. abs(g) ) then
896  c = g/f
897  r = sqrt( c**2 + 1.0 )
898  e(i+1) = f*r
899  s = 1.0/r
900  c = c*s
901  else
902  s = f/g
903  r = sqrt( s**2 + 1.0 )
904  e(i+1) = g*r
905  c = 1.0 / r
906  s = s * c
907  endif
908 c
909  g = d(i+1) - p
910  r = ( d(i) - g ) * s + 2.0 * c * b
911  p = s * r
912  d(i+1) = g + p
913  g = c*r - b
914  14 continue
915 c
916  d(l) = d(l) - p
917  e(l) = g
918  e(m) = 0.0
919  goto 1
920 c
921  endif
922 c
923  15 continue
924 c
925  dmax = 0.0
926  dmin = d(1)
927 c
928  do 40 i = 1 , n
929  dmax = abs( max( d(i) , dmax ) )
930  dmin = abs( min( d(i) , dmin ) )
931  40 continue
932 c
933  return
934  end
935 c-----------------------------------------------------------------------
936  subroutine fdm_h1(z,r,d,mask,mult,nel,kt,rr)
937  include 'SIZE'
938  include 'TOTAL'
939 c
940  common /ctmp0/ w(lx1,ly1,lz1)
941 c
942  include 'FDMH1'
943 c
944 c Overlapping Schwarz, FDM based
945 c
946  real z(lx1,ly1,lz1,1)
947  real r(lx1,ly1,lz1,1)
948  real d(lx1,ly1,lz1,1)
949  real mask(lx1,ly1,lz1,1)
950  real mult(lx1,ly1,lz1,1)
951  real rr(lx1,ly1,lz1,1)
952 c
953  integer kt(lelt,3)
954 c
955  integer icalld
956  save icalld
957  data icalld /0/
958 c
959  n1 = lx1
960  n2 = lx1*lx1
961  n3 = lx1*lx1*lx1
962  ntot = lx1*ly1*lz1*nel
963 c
964  if (ifbhalf) then
965  call col3(rr,r,bhalf,ntot)
966  else
967  call copy(rr,r,ntot)
968 c call col2(rr,mult,ntot)
969  endif
970 c if (nid.eq.0.and.icalld.eq.0) write(6,*) 'In fdm_h1',nel
971  icalld = icalld+1
972 c
973 c
974  do ie=1,nel
975  if (if3d) then
976 c Transfer to wave space:
977  call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n2)
978  do iz=1,n1
979  call mxm(w(1,1,iz),n1,fds(1,kt(ie,2)),n1,z(1,1,iz,ie),n1)
980  enddo
981  call mxm(z(1,1,1,ie),n2,fds(1,kt(ie,3)),n1,w,n1)
982 c
983 c fdsolve:
984 c
985  call col2(w,d(1,1,1,ie),n3)
986 c
987 c Transfer to physical space:
988 c
989  call mxm(w,n2,fdst(1,kt(ie,3)),n1,z(1,1,1,ie),n1)
990  do iz=1,n1
991  call mxm(z(1,1,iz,ie),n1,fdst(1,kt(ie,2)),n1,w(1,1,iz),n1)
992  enddo
993  call mxm(fds(1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n2)
994 c
995  else
996 c Transfer to wave space:
997  call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n1)
998  call mxm(w,n1,fds(1,kt(ie,2)),n1,z(1,1,1,ie),n1)
999 c
1000 c fdsolve:
1001 c
1002  call col2(z(1,1,1,ie),d(1,1,1,ie),n2)
1003 c
1004 c Transfer to physical space:
1005 c
1006  call mxm(z(1,1,1,ie),n1,fdst(1,kt(ie,2)),n1,w,n1)
1007  call mxm(fds(1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n1)
1008 c
1009  endif
1010  enddo
1011 c
1012 c call copy(vx,rr,ntot)
1013 c call copy(vy,z,ntot)
1014 c call prepost(.true.)
1015 c write(6,*) 'quit in fdm'
1016 c call exitt
1017 c
1018  if (ifbhalf) call col2(z,bhalf,ntot)
1019 c
1020 c call col2 (z,mult,ntot)
1021  call dssum(z,lx1,ly1,lz1)
1022  call col2 (z,mask,ntot)
1023 c
1024  return
1025  end
1026 c-----------------------------------------------------------------------
1027  subroutine set_fdm_prec_h1a_gen
1028 c
1029  include 'SIZE'
1030  include 'DXYZ'
1031  include 'INPUT'
1032  include 'MASS'
1033  include 'WZ'
1034 c
1035  include 'FDMH1'
1036 c
1037  COMMON /ctmp0/ w(lx1,lx1),aa(lx1,lx1),bb(lx1,lx1)
1038 c
1039  integer left,right
1040 c
1041 c Set up generic operators for fdm applied to H1 operator (Helmholtz)
1042 c
1043 c 3 cases: E (or P), "D" or "N" for E-E bc, Dirichlet, or Neuamann.
1044 c
1045 c Since there are 2 endpoints, there are a total of 9 types.
1046 c
1047 c
1048  n = lx1
1049  n2 = lx1*lx1
1050 c
1051  delta = abs( zgm1(2,1) - zgm1(1,1) )
1052  bbh = 0.5*delta
1053  aah = 1./delta
1054 c
1055  l = 0
1056  do right = 1,3
1057  do left = 1,3
1058  l = l+1
1059 c
1060  call rzero(bb,n2)
1061  do i=1,lx1
1062  bb(i,i) = wxm1(i)
1063  enddo
1064 c
1065 c A = D^T B D
1066 c
1067  call mxm(bb,n,dxm1 ,n,w,n)
1068  call mxm(dxtm1,n,w,n,aa,n)
1069  if (left.eq.1) then
1070 c Internal
1071  bb(1,1) = bb(1,1) + bbh
1072  aa(1,1) = aa(1,1) + aah
1073  elseif (left.eq.2) then
1074 c Dirichlet
1075  bb(1,1) = 1.
1076  do i=1,n
1077  aa(i,1) = 0.
1078  aa(1,i) = 0.
1079  enddo
1080  aa(1,1) = 1.
1081  endif
1082 c
1083  if (right.eq.1) then
1084 c Internal
1085  bb(n,n) = bb(n,n) + bbh
1086  aa(n,n) = aa(n,n) + aah
1087  elseif (right.eq.2) then
1088 c Dirichlet
1089  bb(n,n) = 1.
1090  do i=1,n
1091  aa(i,n) = 0.
1092  aa(n,i) = 0.
1093  enddo
1094  aa(n,n) = 1.
1095  endif
1096 c
1097 c Scale out mass matrix, so we can precondition w/ binvhf.
1098 c
1099 c ifbhalf = .true.
1100 
1101  ifbhalf = .false.
1102  if (ifbhalf) call rescale_abhalf (aa,bb,w,n)
1103 c
1104 c Now, compute eigenvectors/eigenvalues
1105 c
1106  call generalev(aa,bb,dd(1,l),n,w)
1107  call copy(fds(1,l),aa,n*n)
1108  call transpose(fdst(1,l),n,fds(1,l),n)
1109 c
1110  enddo
1111  enddo
1112  ntot = lx1*ly1*lz1*nelv
1113  if (ifbhalf) call copy (bhalf,binvm1,ntot)
1114  if (ifbhalf) call vsqrt(bhalf,ntot)
1115 c
1116  return
1117  end
1118 c-----------------------------------------------------------------------
1119  subroutine set_fdm_prec_h1a_els
1120 c
1121  include 'SIZE'
1122  include 'DXYZ'
1123  include 'FDMH1'
1124  include 'GEOM'
1125  include 'INPUT'
1126  include 'SOLN'
1127  include 'TOPOL'
1128  include 'WZ'
1129 c
1130  COMMON /ctmp0/ w(lx1,lx1),aa(lx1,lx1),bb(lx1,lx1)
1131  $ , mask(lx1,ly1,lz1,lelt)
1132  real mask
1133  character*3 cb
1134 c
1135 c
1136 c Set up element specific information
1137 c
1138 c 3 cases: E (or P), "D" or "N" for E-E bc, Dirichlet, or Neuamann.
1139 c
1140 c Since there are 2 endpoints, there are a total of 9 types.
1141 c
1142 c
1143  ntot = lx1*ly1*lz1*nelt
1144  kf0 = 1
1145  kf1 = 0
1146  if (ifheat) kf0 = 0
1147  if (ifflow) kf1 = ldim
1148  if (ifsplit) kf1 = ldim+1
1149  do kfld=kf0,kf1
1150  ifld = 1
1151  if (kfld.eq.0) ifld = 2
1152 c
1153  if (kfld.eq.0) call copy(mask, tmask,ntot)
1154  if (kfld.eq.1) call copy(mask,v1mask,ntot)
1155  if (kfld.eq.2) call copy(mask,v2mask,ntot)
1156  if (kfld.eq.3) call copy(mask,v3mask,ntot)
1157  if (kfld.eq.ldim+1) call copy(mask, pmask,ntot)
1158 c
1159  do ie=1,nelv
1160  do ifacedim = 1,ldim
1161 c
1162 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1163 c Mask pointers
1164 c
1165  ii = 2
1166  jj = 2
1167  kk = 2
1168 c
1169  if (ifacedim.eq.1) ii = 1
1170  if (ifacedim.eq.2) jj = 1
1171  if (ifacedim.eq.3) kk = 1
1172  k1 = ii+lx1*(jj-1)
1173  if (if3d) k1 = ii+lx1*(jj-1) + lx1*lx1*(kk-1)
1174 c
1175  if (ifacedim.eq.1) ii = lx1
1176  if (ifacedim.eq.2) jj = lx1
1177  if (ifacedim.eq.3) kk = lx1
1178  k2 = ii+lx1*(jj-1)
1179  if (if3d) k2 = ii+lx1*(jj-1) + lx1*lx1*(kk-1)
1180 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1181 c
1182  iface = 2*ifacedim-1
1183  jface = iface+1
1184 c
1185 c Convert to preproc :(
1186  iface = eface(iface)
1187  jface = eface(jface)
1188 c
1189 c "left" bc
1190 c
1191  cb = cbc(iface,ie,ifld)
1192  if (cb.eq.'E '.or.cb.eq.'P '.or.cb.eq.'p ') then
1193 c Internal
1194  ic1 = 1
1195  elseif (mask(k1,1,1,ie).eq.0) then
1196 c Dirichlet
1197  ic1 = 2
1198  else
1199 c Neumann
1200  ic1 = 3
1201  endif
1202 c write(6,*) ie,iface,'cbl: ',cb,ic1,k1,mask(k1,1,1,ie)
1203 c
1204 c "right" bc
1205 c
1206  cb = cbc(jface,ie,ifld)
1207  if (cb.eq.'E '.or.cb.eq.'P '.or.cb.eq.'p ') then
1208 c Internal
1209  jc1 = 1
1210  elseif (mask(k2,1,1,ie).eq.0) then
1211 c Dirichlet
1212  jc1 = 2
1213  else
1214 c Neumann
1215  jc1 = 3
1216  endif
1217 c write(6,*) ie,jface,'cbr: ',cb,jc1,k2,mask(k2,1,1,ie)
1218 c
1219  ijc = ic1 + 3*(jc1-1)
1220  ktype(ie,ifacedim,kfld) = ijc
1221 c
1222  enddo
1223  enddo
1224  enddo
1225 c
1226 c Boundary condition issues resolved... now resolve length scales
1227 c
1228 c
1229 c
1230  do ie = 1,nelt
1231  do idim=1,ldim
1232  k1 = 1
1233  k2 = lz1
1234  if (idim.eq.3.or.ldim.eq.2) k2=1
1235  j1 = 1
1236  j2 = ly1
1237  if (idim.eq.2) j2=1
1238  i1 = 1
1239  i2 = lx1
1240  if (idim.eq.1) i2=1
1241 c
1242 c l -- face 1, l+jump = face 2
1243 c
1244  jump = (lx1-1)*lx1**(idim-1)
1245  l = 0
1246  dlm = 0
1247  wgt = 0
1248  do k=k1,k2
1249  do j=j1,j2
1250  do i=i1,i2
1251  l = l+1
1252  dl2 = (xm1(i+jump,j,k,ie)-xm1(i,j,k,ie))**2
1253  $ + (ym1(i+jump,j,k,ie)-ym1(i,j,k,ie))**2
1254  $ + (zm1(i+jump,j,k,ie)-zm1(i,j,k,ie))**2
1255  dlm = dlm + dl2*wxm1(i)*wxm1(j)*wxm1(k)
1256  wgt = wgt + wxm1(i)*wxm1(j)*wxm1(k)
1257 c
1258  enddo
1259  enddo
1260  enddo
1261 c
1262  dlm = sqrt(dlm/wgt)
1263  elsize(idim,ie) = dlm/2.
1264 c
1265  enddo
1266 c write(6,1) ie,' elsize:',(elsize(k,ie),k=1,ldim)
1267  enddo
1268  1 format(i8,a8,1p3e15.4)
1269 c
1270  return
1271  end
1272 c-----------------------------------------------------------------------
1273  subroutine set_fdm_prec_h1b(d,h1,h2,nel)
1274  include 'SIZE'
1275  include 'FDMH1'
1276  include 'INPUT'
1277  include 'GEOM'
1278  real d (lx1,ly1,lz1,1)
1279  real h1(lx1,ly1,lz1,1)
1280  real h2(lx1,ly1,lz1,1)
1281 c
1282 c Set up diagonal for FDM for each spectral element
1283 c
1284  nxyz = lx1*ly1*lz1
1285  if (if3d) then
1286  do ie=1,nel
1287  h1b = vlsum(h1(1,1,1,ie),nxyz)/nxyz
1288  h2b = vlsum(h2(1,1,1,ie),nxyz)/nxyz
1289  k1 = ktype(ie,1,kfldfdm)
1290  k2 = ktype(ie,2,kfldfdm)
1291  k3 = ktype(ie,3,kfldfdm)
1292  vol = elsize(1,ie)*elsize(2,ie)*elsize(3,ie)
1293  vl1 = elsize(2,ie)*elsize(3,ie)/elsize(1,ie)
1294  vl2 = elsize(1,ie)*elsize(3,ie)/elsize(2,ie)
1295  vl3 = elsize(1,ie)*elsize(2,ie)/elsize(3,ie)
1296  do i3=1,lz1
1297  do i2=1,ly1
1298  do i1=1,lx1
1299  den = h1b*(vl1*dd(i1,k1) + vl2*dd(i2,k2) + vl3*dd(i3,k3))
1300  $ + h2b*vol
1301  if (ifbhalf) den = den/vol
1302  if (den.ne.0) then
1303  d(i1,i2,i3,ie) = 1./den
1304  else
1305  d(i1,i2,i3,ie) = 0.
1306 c
1307 c write(6,3) 'd=0:'
1308 c $ ,h1(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2),dd(i3,k3)
1309 c $ ,i1,i2,i3,ie,kfldfdm,k1,k2,k3
1310  3 format(a4,1p4e12.4,8i8)
1311 c
1312  endif
1313  enddo
1314  enddo
1315  enddo
1316  enddo
1317  else
1318  do ie=1,nel
1319  if (ifaxis) then
1320  h1b = vlsc2(h1(1,1,1,ie),ym1(1,1,1,ie),nxyz)/nxyz
1321  h2b = vlsc2(h2(1,1,1,ie),ym1(1,1,1,ie),nxyz)/nxyz
1322  else
1323  h1b = vlsum(h1(1,1,1,ie),nxyz)/nxyz
1324  h2b = vlsum(h2(1,1,1,ie),nxyz)/nxyz
1325  endif
1326  k1 = ktype(ie,1,kfldfdm)
1327  k2 = ktype(ie,2,kfldfdm)
1328  vol = elsize(1,ie)*elsize(2,ie)
1329  vl1 = elsize(2,ie)/elsize(1,ie)
1330  vl2 = elsize(1,ie)/elsize(2,ie)
1331  i3=1
1332  do i2=1,ly1
1333  do i1=1,lx1
1334  den = h1b*( vl1*dd(i1,k1) + vl2*dd(i2,k2) )
1335  $ + h2b*vol
1336  if (ifbhalf) den = den/vol
1337  if (den.ne.0) then
1338  d(i1,i2,i3,ie) = 1./den
1339 c write(6,3) 'dn0:'
1340 c $ ,d(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2)
1341 c $ ,i1,i2,i3,ie,kfldfdm,k1,k2
1342  else
1343  d(i1,i2,i3,ie) = 0.
1344 c write(6,3) 'd=0:'
1345 c $ ,h1(i1,i2,i3,ie),dd(i1,k1),dd(i2,k2)
1346 c $ ,i1,i2,i3,ie,kfldfdm,k1,k2
1347  2 format(a4,1p3e12.4,8i8)
1348  endif
1349 c write(6,1) ie,i1,i2,k1,k2,'d:',d(i1,i2,i3,ie),vol,vl1,vl2
1350 c 1 format(5i3,2x,a2,1p4e12.4)
1351  enddo
1352  enddo
1353  enddo
1354  endif
1355 c
1356  return
1357  end
1358 c-----------------------------------------------------------------------
1359  subroutine set_fdm_prec_h1a
1360  include 'SIZE'
1361 c
1362  call set_fdm_prec_h1a_gen
1363  call set_fdm_prec_h1a_els
1364 c
1365  return
1366  end
1367 c-----------------------------------------------------------------------
1368  subroutine generalev(a,b,lam,n,w)
1369 c
1370 c Solve the generalized eigenvalue problem A x = lam B x
1371 c
1372 c A -- symm.
1373 c B -- symm., pos. definite
1374 c
1375 c "SIZE" is included here only to deduce WDSIZE, the working
1376 c precision, in bytes, so as to know whether dsygv or ssygv
1377 c should be called.
1378 c
1379  include 'SIZE'
1380  include 'PARALLEL'
1381 c
1382  real a(n,n),b(n,n),lam(n),w(n,n)
1383  real aa(100),bb(100)
1384 c
1385  parameter(lbw=4*lx1*ly1*lz1*lelv)
1386  common /bigw/ bw(lbw)
1387 c
1388  lw = n*n
1389 c write(6,*) 'in generalev, =',info,n,ninf
1390 c
1391 c call outmat2(a,n,n,n,'aa ')
1392 c call outmat2(b,n,n,n,'bb ')
1393 c
1394  call copy(aa,a,100)
1395  call copy(bb,b,100)
1396 c
1397  call dsygv(1,'V','U',n,a,n,b,n,lam,bw,lbw,info)
1398 c
1399 c call outmat2(a,n,n,n,'Aeig')
1400 c call outmat2(lam,1,n,n,'Deig')
1401 c
1402  if (info.ne.0) then
1403 c
1404  if (nid.eq.0) then
1405  call outmat2(aa ,n,n,n,'aa ')
1406  call outmat2(bb ,n,n,n,'bb ')
1407  call outmat2(a ,n,n,n,'Aeig')
1408  call outmat2(lam,1,n,n,'Deig')
1409  endif
1410 c
1411  ninf = n-info
1412  write(6,*) 'Error in generalev, info=',info,n,ninf
1413  call exitt
1414  endif
1415 c
1416  return
1417  end
1418 c-----------------------------------------------------------------------
1419  subroutine outmat2(a,m,n,k,name)
1420  include 'SIZE'
1421  real a(m,n)
1422  character*4 name
1423 c
1424  n2 = min(n,8)
1425  write(6,2) nid,name,m,n,k
1426  do i=1,m
1427  write(6,1) nid,name,(a(i,j),j=1,n2)
1428  enddo
1429 c 1 format(i3,1x,a4,16f6.2)
1430  1 format(i3,1x,a4,1p8e14.5)
1431  2 format(/,'Matrix: ',i3,1x,a4,3i8)
1432  return
1433  end
1434 c-----------------------------------------------------------------------
1435  subroutine rescale_abhalf (a,b,w,n)
1436  real a(n,n),b(n,n),w(n)
1437 c
1438 c -1/2 -1/2
1439 c Set A = B A B
1440 c
1441 c
1442 c NOTE: B is *diagonal*
1443 c
1444 c
1445  do i=1,n
1446  w(i) = 1./sqrt(b(i,i))
1447  enddo
1448 c
1449  do j=1,n
1450  do i=1,n
1451  a(i,j) = a(i,j)*w(i)*w(j)
1452  enddo
1453  enddo
1454 c
1455 c duh.... don't forget to change B ... duh...
1456 c
1457  call ident(b,n)
1458 c
1459  return
1460  end
1461 c-----------------------------------------------------------------------
1462  subroutine hmholtz_dg(name,u,rhs,h1,h2,mask,tol,maxit)
1463  include 'SIZE'
1464  include 'CTIMER'
1465  include 'INPUT'
1466  include 'MASS'
1467  include 'SOLN'
1468  include 'TSTEP'
1469 C
1470  character name*4
1471  real u (lx1,ly1,lz1,1)
1472  real rhs (lx1,ly1,lz1,1)
1473  real h1 (lx1,ly1,lz1,1)
1474  real h2 (lx1,ly1,lz1,1)
1475  real mask (lx1,ly1,lz1,1)
1476 
1477  icalld=icalld+1
1478  nhmhz=icalld
1479  etime1=dnekclock()
1480 
1481  if (ifield.eq.2) then
1482  call cggo_dg (u,rhs,h1,h2,bintm1,mask,name,tol,maxit)
1483  else
1484  call cggo_dg (u,rhs,h1,h2,binvm1,mask,name,tol,maxit)
1485  endif
1486 
1487  thmhz=thmhz+(dnekclock()-etime1)
1488 
1489  return
1490  end
1491 C
1492 c=======================================================================
1493  subroutine cggo_dg(x,f,h1,h2,binv,mask,name,tin,maxit)
1494 C-------------------------------------------------------------------------
1495 C
1496 C Solve the Helmholtz equation, H*U = RHS,
1497 C using preconditioned conjugate gradient iteration.
1498 C Preconditioner: diag(H).
1499 C
1500 C------------------------------------------------------------------------
1501  include 'SIZE'
1502  include 'TOTAL'
1503 
1504 
1505  real x(1),f(1),h1(1),h2(1),binv(1),mask(1)
1506  parameter(lg=lx1*ly1*lz1*lelt)
1507  common /scrcg/ d(lg) , scalar(2)
1508  common /scrmg/ r(lg) , w(lg) , p(lg) , z(lg)
1509 
1510  parameter(maxcg=900)
1511  common /tdarray/ diagt(maxcg),upper(maxcg)
1512  common /iterhm/ niterhm
1513  character*4 name
1514 
1515  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1516  logical ifdfrm, iffast, ifh2, ifsolv
1517 
1518  common /cprint/ ifprint, ifhzpc
1519  logical ifprint, ifhzpc
1520 
1521  logical ifmcor
1522 
1523 
1524 c ** zero out stuff for Lanczos eigenvalue estimator
1525  call rzero(diagt,maxcg)
1526  call rzero(upper,maxcg)
1527 
1528 c Initialization
1529 
1530  nxyz = lx1*ly1*lz1
1531  nel = nelv
1532  vol = volvm1
1533  if (ifield.eq.2) nel=nelt
1534  if (ifield.eq.2) vol=voltm1
1535  n = nel*nxyz
1536 
1537  tol=tin
1538  if (param(22).ne.0) tol=abs(param(22))
1539  niter = min(maxit,maxcg)
1540 
1541  imsh = ifield
1542  call setprec_dg(d,h1,h2,imsh,1) ! diag preconditioner
1543 c call invers2 (d,bm1,n) ! diag preconditioner
1544 
1545  call copy (r,f,n)
1546  call rzero(x,n)
1547  call rzero(p,n)
1548 
1549 c Check for non-trivial null-space
1550 
1551  ifmcor = .false.
1552  h2max = glmax(h2 ,n)
1553  skmin = glmin(mask,n)
1554  if (skmin.gt.0.and.h2max.eq.0) ifmcor = .true.
1555 
1556  if (ifmcor) then
1557  rmean = glsum(r,n)
1558  call cadd(r,rmean,n)
1559  endif
1560 
1561  krylov = 0
1562  rtz1=1.0
1563  niterhm = 0
1564  do 1000 iter=1,niter
1565 
1566 c call copy(z,r,n) ! No preconditioner
1567  call col3(z,r,d,n) ! Jacobi Preconditioner
1568 
1569  rtz2=rtz1
1570  scalar(1)=vlsc2(z,r,n)
1571  scalar(2)=vlsc3(r,r,binv,n)
1572  call gop(scalar,w,'+ ',2)
1573  rtz1=scalar(1)
1574  rbn2=sqrt(scalar(2)/vol)
1575  if (iter.eq.1) rbn0 = rbn2
1576  if (param(22).lt.0) tol=abs(param(22))*rbn0
1577 
1578  if (ifprint.and.nid.eq.0.and.param(74).ne.0) then
1579  write(6,3002) istep,iter,name,ifmcor,rbn2,tol,h1(1),h2(1)
1580  endif
1581 
1582  if (rbn2.le.tol) then
1583  niter = iter-1
1584  if(nid.eq.0.and.((.not.ifhzpc).or.ifprint))
1585  $ write(6,3000) istep,name,niter,rbn2,rbn0,tol
1586  go to 9999
1587  endif
1588 
1589  beta = rtz1/rtz2
1590  if (iter.eq.1) beta=0.0
1591  call add2s1 (p,z,beta,n)
1592  call hxdg (w,p,h1,h2)
1593 
1594  rho0 = rho
1595  rho = glsc2(w,p,n)
1596  alpha=rtz1/rho
1597  alphm=-alpha
1598  call add2s2(x,p ,alpha,n)
1599  call add2s2(r,w ,alphm,n)
1600 
1601 c Generate tridiagonal matrix for Lanczos scheme
1602  if (iter.eq.1) then
1603  krylov = krylov+1
1604  diagt(iter) = rho/rtz1
1605  elseif (iter.le.maxcg) then
1606  krylov = krylov+1
1607  diagt(iter) = (beta**2 * rho0 + rho ) / rtz1
1608  upper(iter-1) = -beta * rho0 / sqrt(rtz2 * rtz1)
1609  endif
1610  1000 continue
1611  niter = iter-1
1612 c
1613  if (nid.eq.0) write (6,3001) istep,niter,name,rbn2,rbn0,tol
1614  3000 format(i9,4x,'hmh dg ',a4,': ',i6,1p6e13.4)
1615  3001 format(2i6,' **ERROR**: Failed in hmh_dg: ',a4,1p6e13.4)
1616  3002 format(i3,i6,' HMH dg ',a4,1x,l4,':',1p6e13.4)
1617  9999 continue
1618  niterhm = niter
1619  ifsolv = .false.
1620 c
1621 c
1622 c Call eigenvalue routine for Lanczos scheme:
1623 c two work arrays are req'd if you want to save "diag & upper"
1624 c
1625 c if (iter.ge.3) then
1626 c niter = iter-1
1627 c call calc (diagt,upper,w,z,krylov,dmax,dmin)
1628 c cond = dmax/dmin
1629 c if (nid.eq.0) write(6,6) istep,cond,dmin,dmax,' lambda'
1630 c endif
1631 c 6 format(i9,1p3e12.4,4x,a7)
1632 c
1633 c if (n.gt.0) write(6,*) 'quit in cggo'
1634 c if (n.gt.0) call exitt
1635 c call exitt
1636  return
1637  end
1638 c-----------------------------------------------------------------------
1639  subroutine outmax(a,m,n,name6,ie)
1640  real a(m,n)
1641  character*6 name6
1642 c
1643  n18 = min(n,18)
1644  write(6,*)
1645  write(6,*) ie,' matrix: ',name6,m,n
1646  do i=1,m
1647  write(6,6) ie,name6,(a(i,j),j=1,n18)
1648  enddo
1649  6 format(i3,1x,a6,18f7.2)
1650  write(6,*)
1651  return
1652  end
1653 c-----------------------------------------------------------------------
1654  subroutine outmat4(a,l,m,n,nel,name6,ie)
1655  real a(l,m,n,nel)
1656  character*6 name6
1657 c
1658  n18 = min(n,18)
1659  write(6,*)
1660  write(6,*) ie,' matrix: ',name6,m,n
1661  do ll=1,l
1662  do k=1,nel
1663  write(6,*) ie,' matrix: ',name6,ll,k
1664  do j=1,4
1665  write(6,6) ie,name6,(a(ll,i,j,k),i=1,m)
1666  enddo
1667  enddo
1668  enddo
1669  6 format(i3,1x,a6,18f7.2)
1670  write(6,*)
1671  return
1672  end
1673 c-----------------------------------------------------------------------
1674  subroutine ioutmat4(a,l,m,n,nel,name6,ie)
1675  integer a(l,m,n,nel)
1676  character*6 name6
1677 c
1678  n18 = min(n,18)
1679  write(6,*)
1680  write(6,*) ie,' matrix: ',name6,m,n
1681  do ll=1,l
1682  do k=1,nel
1683  write(6,*) ie,' matrix: ',name6,ll,k
1684  do j=1,n
1685  write(6,6) ie,name6,(a(ll,i,j,k),i=1,m)
1686  enddo
1687  enddo
1688  enddo
1689  6 format(i3,1x,a6,18i7)
1690  write(6,*)
1691  return
1692  end
1693 c-----------------------------------------------------------------------
1694  subroutine ioutfld(a,m,n,nel,name6,ie)
1695  integer a(m,n,nel)
1696  character*6 name6
1697 
1698  n18 = min(n,18)
1699  write(6,*)
1700  write(6,*) ie,' matrix: ',name6,m,n
1701  do j=1,n
1702  if (m.eq.3) write(6,3) ie,name6,((a(i,j,k),i=1,m),k=1,2)
1703  if (m.eq.4) write(6,4) ie,name6,((a(i,j,k),i=1,m),k=1,2)
1704  if (m.eq.5) write(6,5) ie,name6,((a(i,j,k),i=1,m),k=1,2)
1705  if (m.eq.6) write(6,6) ie,name6,((a(i,j,k),i=1,m),k=1,2)
1706  enddo
1707  3 format(i3,1x,a6,2(3i7,2x))
1708  4 format(i3,1x,a6,2(4i7,2x))
1709  5 format(i3,1x,a6,2(5i7,2x))
1710  6 format(i3,1x,a6,2(6i7,2x))
1711  write(6,*)
1712  return
1713  end
1714 c-----------------------------------------------------------------------
1715  subroutine gradr(ur,us,ut,u,Dr,Dst,Dtt,nr,ns,nt,if3d)
1716 c
1717 c Output: ur,us,ut Input: u
1718 c
1719  real ur(nr,ns,nt),us(nr,ns,nt),ut(nr,ns,nt)
1720  real u (nr,ns,nt)
1721  real Dr(nr,nr),Dst(ns,ns),Dtt(nt,nt)
1722 c
1723  logical if3d
1724 c
1725  nst = ns*nt
1726  nrs = nr*ns
1727 c
1728  if (if3d) then
1729  call mxm(dr,nr,u,nr,ur,nst)
1730  do k=1,nt
1731  call mxm(u(1,1,k),nr,dst,ns,us(1,1,k),nt)
1732  enddo
1733  call mxm(u,nrs,dtt,nt,ut,nt)
1734  else
1735  call mxm(dr,nr,u ,nr,ur,ns)
1736  call mxm(u ,nr,dst,ns,us,ns)
1737  endif
1738 c
1739  return
1740  end
1741 c-----------------------------------------------------------------------
1742  subroutine gradrta(u,ur,us,ut,Drt,Ds,Dt,nr,ns,nt,if3d)
1743 c
1744 c T T T
1745 c Output: u = u + D ur + D us + D ut Input: ur,us,ut
1746 c r s t
1747 c
1748  real u (nr,ns,nt)
1749  real ur(nr,ns,nt),us(nr,ns,nt),ut(nr,ns,nt)
1750  real Dr(nr,nr),Dst(ns,ns),Dtt(nt,nt)
1751 c
1752  logical if3d
1753 c
1754  nst = ns*nt
1755  nrs = nr*ns
1756 c
1757  if (if3d) then
1758  call mxma(drt,nr,ur,nr,u,nst)
1759  do k=1,nt
1760  call mxma(us(1,1,k),nr,ds,ns,u(1,1,k),nt)
1761  enddo
1762  call mxma(ut,nrs,dt,nt,u,nt)
1763  else
1764  call mxma(drt,nr,ur,nr,u,ns)
1765  call mxma(us ,nr,ds,ns,u,ns)
1766  endif
1767 c
1768  return
1769  end
1770 c-----------------------------------------------------------------------
1771  subroutine face_diff (u,d,gsh_loc,w) ! difference: e_f - e'_f
1772 
1773  include 'SIZE'
1774  include 'TOPOL'
1775  include 'PARALLEL'
1776 
1777  integer d,gsh_loc
1778  real u(lx1*lz1*2*ldim*lelt,2)
1779  real w(lx1*lz1*2*ldim*lelt,2)
1780 
1781  n = 2*ldim*lx1*lz1*nelt
1782 
1783  do j=1,d
1784  do i=1,n
1785  w(i,j) = u(i,j)
1786  enddo
1787  call fgslib_gs_op (gsh_loc,w(1,j),1,1,0) ! 1 ==> +
1788 
1789  do i=1,n
1790  u(i,j) = 2*u(i,j)-w(i,j)
1791  enddo
1792 
1793  enddo
1794 
1795  return
1796  end
1797 c-----------------------------------------------------------------------
1798  subroutine setprec_dg (d,h1,h2,imsh,isd)
1799 C-------------------------------------------------------------------
1800 C
1801 C Generate diagonal preconditioner for the DG Helmholtz operator.
1802 C
1803 C-------------------------------------------------------------------
1804  include 'SIZE'
1805  include 'TOTAL'
1806  real d(lx1,ly1,lz1,1)
1807  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1808  logical ifdfrm, iffast, ifh2, ifsolv
1809  real h1(lx1*ly1*lz1,1), h2(lx1*ly1*lz1,1)
1810  real ysm1(ly1)
1811  integer e,f,pf
1812 
1813  nel=nelt
1814  if (imsh.eq.1) nel=nelv
1815 
1816  call dsset(lx1,ly1,lz1)
1817 
1818  n = nel*lx1*ly1*lz1
1819  nxyz = lx1*ly1*lz1
1820  nface = 2*ldim
1821 
1822  do 1000 e=1,nel
1823 
1824  call rzero(d(1,1,1,e),nxyz)
1825 
1826  if (ldim.eq.3) then
1827 
1828  do 320 iz=1,lz1
1829  do 320 iy=1,ly1
1830  do 320 ix=1,lx1
1831  do 320 iq=1,lx1
1832  d(ix,iy,iz,e) = d(ix,iy,iz,e)
1833  $ + g1m1(iq,iy,iz,e) * dxm1(iq,ix)**2
1834  $ + g2m1(ix,iq,iz,e) * dxm1(iq,iy)**2
1835  $ + g3m1(ix,iy,iq,e) * dxm1(iq,iz)**2
1836  320 continue
1837 c
1838 c Add cross terms if element is deformed.
1839 c
1840  if (ifdfrm(e)) then
1841 
1842  do i2=1,ly1,ly1-1
1843  do i1=1,lx1,lx1-1
1844  d(1,i1,i2,e) = d(1,i1,i2,e)
1845  $ + g4m1(1,i1,i2,e) * dxtm1(1,1)*dytm1(i1,i1)
1846  $ + g5m1(1,i1,i2,e) * dxtm1(1,1)*dztm1(i2,i2)
1847  d(lx1,i1,i2,e) = d(lx1,i1,i2,e)
1848  $ + g4m1(lx1,i1,i2,e) * dxtm1(lx1,lx1)*dytm1(i1,i1)
1849  $ + g5m1(lx1,i1,i2,e) * dxtm1(lx1,lx1)*dztm1(i2,i2)
1850  d(i1,1,i2,e) = d(i1,1,i2,e)
1851  $ + g4m1(i1,1,i2,e) * dytm1(1,1)*dxtm1(i1,i1)
1852  $ + g6m1(i1,1,i2,e) * dytm1(1,1)*dztm1(i2,i2)
1853  d(i1,ly1,i2,e) = d(i1,ly1,i2,e)
1854  $ + g4m1(i1,ly1,i2,e) * dytm1(ly1,ly1)*dxtm1(i1,i1)
1855  $ + g6m1(i1,ly1,i2,e) * dytm1(ly1,ly1)*dztm1(i2,i2)
1856  d(i1,i2,1,e) = d(i1,i2,1,e)
1857  $ + g5m1(i1,i2,1,e) * dztm1(1,1)*dxtm1(i1,i1)
1858  $ + g6m1(i1,i2,1,e) * dztm1(1,1)*dytm1(i2,i2)
1859  d(i1,i2,lz1,e) = d(i1,i2,lz1,e)
1860  $ + g5m1(i1,i2,lz1,e) * dztm1(lz1,lz1)*dxtm1(i1,i1)
1861  $ + g6m1(i1,i2,lz1,e) * dztm1(lz1,lz1)*dytm1(i2,i2)
1862 
1863  enddo
1864  enddo
1865  endif
1866 
1867  else ! 2d
1868 
1869  iz=1
1870  if (ifaxis) call setaxdy ( ifrzer(e) )
1871 
1872  do 220 iy=1,ly1
1873  do 220 ix=1,lx1
1874  do 220 iq=1,lx1
1875  d(ix,iy,iz,e) = d(ix,iy,iz,e)
1876  $ + g1m1(iq,iy,iz,e) * dxm1(iq,ix)**2
1877  $ + g2m1(ix,iq,iz,e) * dxm1(iq,iy)**2
1878  220 continue
1879 c
1880 
1881  if (ifdfrm(e)) then
1882 
1883  do i1=1,ly1,ly1-1
1884  d(1,i1,iz,e) = d(1,i1,iz,e)
1885  $ + g4m1(1,i1,iz,e) * dxm1(1,1)*dym1(i1,i1)
1886  d(lx1,i1,iz,e) = d(lx1,i1,iz,e)
1887  $ + g4m1(lx1,i1,iz,e) * dxm1(lx1,lx1)*dym1(i1,i1)
1888  d(i1,1,iz,e) = d(i1,1,iz,e)
1889  $ + g4m1(i1,1,iz,e) * dym1(1,1)*dxm1(i1,i1)
1890  d(i1,ly1,iz,e) = d(i1,ly1,iz,e)
1891  $ + g4m1(i1,ly1,iz,e) * dym1(ly1,ly1)*dxm1(i1,i1)
1892  enddo
1893  endif
1894 
1895  endif
1896 
1897 c Here, we add DG surface terms (11/06/16)
1898 
1899  do f=1,nface
1900  pf = eface1(f)
1901  js1 = skpdat(1,pf)
1902  jf1 = skpdat(2,pf)
1903  jskip1 = skpdat(3,pf)
1904  js2 = skpdat(4,pf)
1905  jf2 = skpdat(5,pf)
1906  jskip2 = skpdat(6,pf)
1907 
1908  i = 0
1909  do j2=js2,jf2,jskip2
1910  do j1=js1,jf1,jskip1
1911  i = i+1
1912  d(j1,j2,1,e) = d(j1,j2,1,e) + etalph(i,f,e)
1913  enddo
1914  enddo
1915  enddo
1916 
1917  i=0
1918  nx=lx1
1919  if (ldim.eq.3) then
1920  do i2=1,ly1
1921  do i1=1,lx1
1922  i=i+1
1923  d( 1,i1,i2,e)=d( 1,i1,i2,e)-2*fw(4,e)*unr(i,4,e)*dxm1( 1, 1)
1924  d(nx,i1,i2,e)=d(nx,i1,i2,e)-2*fw(2,e)*unr(i,2,e)*dxm1(nx,nx)
1925  d(i1, 1,i2,e)=d(i1, 1,i2,e)-2*fw(1,e)*uns(i,1,e)*dym1( 1, 1)
1926  d(i1,nx,i2,e)=d(i1,nx,i2,e)-2*fw(3,e)*uns(i,3,e)*dym1(nx,nx)
1927  d(i1,i2, 1,e)=d(i1,i2, 1,e)-2*fw(5,e)*unt(i,5,e)*dzm1( 1, 1)
1928  d(i1,i2,nx,e)=d(i1,i2,nx,e)-2*fw(6,e)*unt(i,6,e)*dzm1(nx,nx)
1929  enddo
1930  enddo
1931  else ! 2D
1932  do i1=1,lx1
1933  i=i+1
1934  d( 1,i1,1,e)=d( 1,i1,1,e)-2*fw(4,e)*unr(i,4,e)*dxm1( 1, 1)
1935  d(nx,i1,1,e)=d(nx,i1,1,e)-2*fw(2,e)*unr(i,2,e)*dxm1(nx,nx)
1936  d(i1, 1,1,e)=d(i1, 1,1,e)-2*fw(1,e)*uns(i,1,e)*dym1( 1, 1)
1937  d(i1,nx,1,e)=d(i1,nx,1,e)-2*fw(3,e)*uns(i,3,e)*dym1(nx,nx)
1938  enddo
1939  endif
1940 
1941  do i=1,nxyz
1942  d(i,1,1,e)=1./(d(i,1,1,e)*h1(i,e)+h2(i,e)*bm1(i,1,1,e))
1943  enddo
1944 
1945  1000 continue ! element loop
1946 
1947 c If axisymmetric, add a diagonal term in the radial direction (ISD=2)
1948 
1949  if (ifaxis.and.(isd.eq.2)) then
1950  call invcol1 (d,n)
1951  do 1200 e=1,nel
1952 
1953  if (ifrzer(e)) call mxm(ym1(1,1,1,e),lx1,datm1,ly1,ysm1,1)
1954 
1955  k=0
1956  do 1190 j=1,ly1
1957  do 1190 i=1,lx1
1958  k=k+1
1959  if (ym1(i,j,1,e).ne.0.) then
1960  term1 = bm1(i,j,1,e)/ym1(i,j,1,e)**2
1961  if (ifrzer(e)) then
1962  term2 = wxm1(i)*wam1(1)*dam1(1,j)
1963  $ *jacm1(i,1,1,e)/ysm1(i)
1964  else
1965  term2 = 0.
1966  endif
1967  d(i,j,1,e) = d(i,j,1,e)
1968  $ + h1(k,e)*(term1+term2)
1969  endif
1970  1190 continue
1971  1200 continue
1972 
1973  call invcol1 (d,n)
1974 
1975  endif
1976 
1977  return
1978  end
1979 c-----------------------------------------------------------------------
1980  subroutine hxdg_surfa(au,u,h1,h2)
1981 
1982 c Helmholtz matrix-vector product: Au = Au + surface term
1983 
1984  include 'SIZE'
1985  include 'TOTAL'
1986 
1987  parameter(lxyz=lx1*ly1*lz1)
1988 
1989  real au(lx1,ly1,lz1,lelt),u(lx1,ly1,lz1,lelt)
1990  real h1(lx1,ly1,lz1,lelt),h2(1)
1991 
1992  common /ytmp9/ qr(lx1,ly1,lz1),qs(lx1,ly1,lz1),qt(lx1,ly1,lz1)
1993 
1994  integer e,f,pf
1995 
1996 
1997  call dsset(lx1,ly1,lz1)
1998  nface = 2*ldim
1999  n = lx1*ly1*lz1*nelfld(ifield)
2000 
2001  do e=1,nelfld(ifield)
2002  iflag=0
2003  do f=1,nface
2004  if (fw(f,e).gt.0.6) iflag=1
2005  enddo
2006  if (iflag.gt.0) then
2007 
2008  if (ifaxis) call setaxdy(ifrzer(e))
2009 
2010  do i=1,lxyz
2011  qr(i,1,1)=0
2012  qs(i,1,1)=0
2013  qt(i,1,1)=0
2014  enddo
2015 
2016  do f=1,nface
2017  if (fw(f,e).gt.0.6) then
2018  pf = eface1(f)
2019  js1 = skpdat(1,pf)
2020  jf1 = skpdat(2,pf)
2021  jskip1 = skpdat(3,pf)
2022  js2 = skpdat(4,pf)
2023  jf2 = skpdat(5,pf)
2024  jskip2 = skpdat(6,pf)
2025 
2026  fwtbc=1.
2027  if (cbc(f,e,ifield).eq.'O '.or.cbc(f,e,ifield).eq.'I ')
2028  $ fwtbc=0
2029 
2030  i = 0
2031  do j2=js2,jf2,jskip2
2032  do j1=js1,jf1,jskip1
2033  i = i+1
2034  fwt = fwtbc * h1(j1,j2,1,e)*u(j1,j2,1,e)
2035  et1 = etalph(i,f,e)*h1(j1,j2,1,e)*u(j1,j2,1,e)
2036  qr(j1,j2,1) = qr(j1,j2,1)-fwt*unr(i,f,e)
2037  qs(j1,j2,1) = qs(j1,j2,1)-fwt*uns(i,f,e)
2038  qt(j1,j2,1) = qt(j1,j2,1)-fwt*unt(i,f,e)
2039  au(j1,j2,1,e) = au(j1,j2,1,e)+et1*fwtbc
2040  enddo
2041  enddo
2042  endif
2043  enddo
2044 
2045  call gradrta(au(1,1,1,e),qr,qs,qt ! NOTE FIX in gradr()! 3D!
2046  $ ,dxtm1,dym1,dzm1,lx1,ly1,lz1,if3d)
2047  endif
2048  enddo
2049 
2050  return
2051  end
2052 c-----------------------------------------------------------------------
2053  subroutine hxdg (au,u,h1,h2)
2054 
2055 c Helmholtz matrix-vector product: Au = h1*[A]u + h2*[B]u
2056 
2057  include 'SIZE'
2058  include 'TOTAL'
2059 
2060  parameter(lxyz=lx1*ly1*lz1)
2061  real au(lx1,ly1,lz1,1),u(lx1,ly1,lz1,1),h1(lx1,ly1,lz1,1),h2(1)
2062 
2063  common /ctmp0/ w(2*lx1*lz1*2*ldim*lelt)
2064  common /ctmp1/ ur(lx1,ly1,lz1,lelt),us(lx1,ly1,lz1,lelt)
2065  $ ,ut(lx1,ly1,lz1,lelt)
2066  common /ytmp9/ qr(lx1,ly1,lz1),qs(lx1,ly1,lz1),qt(lx1,ly1,lz1)
2067  common /ytmp0/ uf(lx1*lz1,2*ldim,lelt,2)
2068 
2069  integer e,f,pf
2070 
2071  n = lx1*ly1*lz1*nelfld(ifield)
2072  nface = 2*ldim
2073 
2074  call dsset(lx1,ly1,lz1)
2075 
2076  call col4(au,h2,bm1,u,n) ! au = h2 B u
2077 
2078  do e=1,nelfld(ifield)
2079 
2080  if (ifaxis) call setaxdy(ifrzer(e))
2081 
2082  call gradr(ur(1,1,1,e),us(1,1,1,e),ut(1,1,1,e) ! NOTE FIX in gradr()! 3D!
2083  $ ,u(1,1,1,e),dxm1,dytm1,dztm1,lx1,ly1,lz1,if3d)
2084 
2085  do f=1,nface
2086  pf = eface1(f)
2087  js1 = skpdat(1,pf)
2088  jf1 = skpdat(2,pf)
2089  jskip1 = skpdat(3,pf)
2090  js2 = skpdat(4,pf)
2091  jf2 = skpdat(5,pf)
2092  jskip2 = skpdat(6,pf)
2093 
2094  i = 0
2095  do j2=js2,jf2,jskip2
2096  do j1=js1,jf1,jskip1
2097  i = i+1
2098 c Normally, we'd store this as a 2-vector: uf(2,...)
2099  uf(i,f,e,1) = u(j1,j2,1,e)*h1(j1,j2,1,e)
2100  uf(i,f,e,2) = (unr(i,f,e)*ur(j1,j2,1,e)
2101  $ + uns(i,f,e)*us(j1,j2,1,e)
2102  $ + unt(i,f,e)*ut(j1,j2,1,e))*h1(j1,j2,1,e)
2103  enddo
2104  enddo
2105  enddo
2106  enddo
2107 
2108  call face_diff (uf,2,dg_hndlx,w) ! difference: e_f - e'_f
2109 
2110  do e=1,nelfld(ifield)
2111  if (ifaxis) call setaxdy(ifrzer(e))
2112  do i=1,lxyz
2113  qr(i,1,1) = (g1m1(i,1,1,e)*ur(i,1,1,e)
2114  $ +g4m1(i,1,1,e)*us(i,1,1,e)
2115  $ +g5m1(i,1,1,e)*ut(i,1,1,e))*h1(i,1,1,e)
2116  qs(i,1,1) = (g4m1(i,1,1,e)*ur(i,1,1,e)
2117  $ +g2m1(i,1,1,e)*us(i,1,1,e)
2118  $ +g6m1(i,1,1,e)*ut(i,1,1,e))*h1(i,1,1,e)
2119  qt(i,1,1) = (g5m1(i,1,1,e)*ur(i,1,1,e)
2120  $ +g6m1(i,1,1,e)*us(i,1,1,e)
2121  $ +g3m1(i,1,1,e)*ut(i,1,1,e))*h1(i,1,1,e)
2122  enddo
2123 
2124  do f=1,nface
2125 
2126  fwtbc=1.
2127  if (cbc(f,e,ifield).eq.'O '.or.cbc(f,e,ifield).eq.'I ')
2128  $ fwtbc=0
2129  pf = eface1(f)
2130  js1 = skpdat(1,pf)
2131  jf1 = skpdat(2,pf)
2132  jskip1 = skpdat(3,pf)
2133  js2 = skpdat(4,pf)
2134  jf2 = skpdat(5,pf)
2135  jskip2 = skpdat(6,pf)
2136 
2137 
2138  i = 0
2139  do j2=js2,jf2,jskip2
2140  do j1=js1,jf1,jskip1
2141  i = i+1
2142  fwt = fw(f,e)*fwtbc
2143  qr(j1,j2,1)=qr(j1,j2,1)-fwt*unr(i,f,e)*uf(i,f,e,1)
2144  qs(j1,j2,1)=qs(j1,j2,1)-fwt*uns(i,f,e)*uf(i,f,e,1)
2145  qt(j1,j2,1)=qt(j1,j2,1)-fwt*unt(i,f,e)*uf(i,f,e,1)
2146  au(j1,j2,1,e) = au(j1,j2,1,e)-fwt*uf(i,f,e,2)
2147  $ + etalph(i,f,e)*uf(i,f,e,1)*fwtbc
2148  enddo
2149  enddo
2150 
2151  enddo
2152 
2153  call gradrta(au(1,1,1,e),qr,qs,qt ! NOTE FIX in gradr()! 3D!
2154  $ ,dxtm1,dym1,dzm1,lx1,ly1,lz1,if3d)
2155  enddo
2156 
2157  return
2158  end
2159 c-----------------------------------------------------------------------
2160  subroutine hmh_flex_cg(res,h1,h2,wt,iter)
2161 
2162 c Solve the Helmholtz equation by right-preconditioned
2163 c GMRES iteration.
2164 
2165 
2166  include 'SIZE'
2167  include 'TOTAL'
2168  include 'FDMH1'
2169  include 'GMRES'
2170  common /ctolpr/ divex
2171  common /cprint/ ifprint
2172  logical ifprint
2173  real res (lx1*ly1*lz1*lelv)
2174  real h1 (lx1,ly1,lz1,lelv)
2175  real h2 (lx1,ly1,lz1,lelv)
2176  real wt (lx1,ly1,lz1,lelv)
2177 
2178  parameter(lt=lx1*ly1*lz1*lelt)
2179  common /scrcg/ r(lt),z(lt),p(lt),w(lt)
2180  common /scrmg/ r1(lt)
2181 
2182  common /cgmres1/ y(lgmres)
2183  common /ctmp0/ wk1(lgmres),wk2(lgmres)
2184  real alpha, l, temp
2185  integer outer
2186 
2187  logical iflag,if_hyb
2188  save iflag,if_hyb
2189 c data iflag,if_hyb /.false. , .true. /
2190  data iflag,if_hyb /.false. , .false. /
2191  real norm_fac
2192  save norm_fac
2193 
2194  real*8 etime1,dnekclock
2195 
2196  n = lx1*ly1*lz1*nelv
2197 
2198  div0 = gamma_gmres(1)*norm_fac
2199 
2200  etime1 = dnekclock()
2201  etime_p = 0.
2202  divex = 0.
2203  maxit = iter
2204  iter = 0
2205 
2206 
2207  call chktcg1(tolps,res,h1,h2,pmask,vmult,1,1)
2208  if (param(21).gt.0.and.tolps.gt.abs(param(21)))
2209  $ tolps = abs(param(21))
2210  if (istep.eq.0) tolps = 1.e-4
2211  tolpss = tolps
2212 
2213  iconv = 0
2214  call copy (r,res,n) ! Residual
2215  call rzero(r1 ,n) ! Lagged residual for flexible CG
2216  call rzero(p,n) ! Search direction
2217  call rzero(res,n) ! Solution vector
2218  rho1 = 1
2219  ! ______
2220  div0 = sqrt(glsc3(r,wt,r,n)/volvm1) ! gamma = \/ (r,r)
2221 
2222  if (param(21).lt.0) tolpss=abs(param(21))*div0
2223 
2224 
2225  do k=1,maxit
2226 
2227  if (param(40).ge.0 .and. param(40).le.2) then
2228  call h1mg_solve(z,r,if_hyb) ! z = M w
2229  else if (param(40).eq.3) then
2230  call fem_amg_solve(z,r)
2231  endif
2232 
2233  call sub2(r1,r,n)
2234  rho0 = rho1
2235  rho1 = glsc3(z,wt,r,n) ! Inner product weighted by multiplicity
2236  rho2 = -glsc3(z,wt,r1,n) ! Inner product weighted by multiplicity
2237  beta = rho2/rho0 ! Flexible GMRES
2238 
2239  call copy(r1,r,n) ! Save prior residual
2240  call add2s1(p,z,beta,n)
2241 
2242  call ax(w,p,h1,h2,n) ! w = A p
2243  den = glsc3(w,wt,p,n)
2244  alpha = rho1/den
2245  rnorm = 0.0
2246  do i = 1,n
2247  res(i) = res(i) + alpha*p(i)
2248  r(i) = r(i) - alpha*w(i)
2249  rnorm = rnorm + r(i)*r(i)*wt(i,1,1,1)
2250  enddo
2251  call gop(rnorm,temp,'+ ',1)
2252 
2253 c rnorm = sqrt(glsc3(r,wt,r,n)/volvm1) ! gamma = \/ (r,r)
2254  rnorm = sqrt(rnorm/volvm1) ! gamma = \/ (r,r)
2255  ratio = rnorm/div0
2256 
2257  iter=iter+1
2258  if (ifprint.and.nio.eq.0)
2259  $ write (6,66) iter,tolpss,rnorm,div0,ratio,istep
2260  66 format(i5,1p4e12.5,i8,' Divergence')
2261 
2262  if (rnorm .lt. tolpss) goto 900 !converged
2263 
2264  enddo
2265 
2266  900 iconv = 1
2267  divex = rnorm
2268  call ortho (res) ! Orthogonalize wrt null space, if present
2269  etime1 = dnekclock()-etime1
2270  if (nio.eq.0) write(6,9999) istep,iter,divex,div0,tolpss,etime_p,
2271  & etime1,if_hyb
2272  9999 format(4x,i7,' PRES cgflex',4x,i5,1p5e13.4,1x,l4)
2273 
2274  if (outer.le.2) if_hyb = .false.
2275 
2276  return
2277  end
2278 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
ident
subroutine ident(a, n)
Definition: calcz.f:147
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
dsygv
subroutine dsygv(ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO)
Definition: dsygv.f:3
fem_amg_solve
#define fem_amg_solve
Definition: fem_amg_preco.c:9
h1_overlap_2
subroutine h1_overlap_2(u, v, mask)
Definition: gmres.f:572
hmh_gmres
subroutine hmh_gmres(res, h1, h2, wt, iter)
Definition: gmres.f:305
ax
subroutine ax(w, x, h1, h2, n)
Definition: gmres.f:285
fdm_h1
subroutine fdm_h1(z, r, d, mask, mult, nel, kt, rr)
Definition: hmholtz.f:937
outmat2
subroutine outmat2(a, m, n, k, name)
Definition: hmholtz.f:1420
setfast
subroutine setfast(helm1, helm2, imesh)
Definition: hmholtz.f:263
chktcg1
subroutine chktcg1(tol, res, h1, h2, mask, mult, imesh, isd)
Definition: hmholtz.f:528
sfastax
subroutine sfastax
Definition: hmholtz.f:311
outmax
subroutine outmax(a, m, n, name6, ie)
Definition: hmholtz.f:1640
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
generalev
subroutine generalev(a, b, lam, n, w)
Definition: hmholtz.f:1369
gradr
subroutine gradr(ur, us, ut, u, Dr, Dst, Dtt, nr, ns, nt, if3d)
Definition: hmholtz.f:1716
set_fdm_prec_h1a_els
subroutine set_fdm_prec_h1a_els
Definition: hmholtz.f:1120
hxdg_surfa
subroutine hxdg_surfa(au, u, h1, h2)
Definition: hmholtz.f:1981
hxdg
subroutine hxdg(au, u, h1, h2)
Definition: hmholtz.f:2054
gradrta
subroutine gradrta(u, ur, us, ut, Drt, Ds, Dt, nr, ns, nt, if3d)
Definition: hmholtz.f:1743
cggo
subroutine cggo(x, f, h1, h2, mask, mult, imsh, tin, maxit, isd, binv, name)
Definition: hmholtz.f:612
hmh_flex_cg
subroutine hmh_flex_cg(res, h1, h2, wt, iter)
Definition: hmholtz.f:2161
outmat4
subroutine outmat4(a, l, m, n, nel, name6, ie)
Definition: hmholtz.f:1655
ioutmat4
subroutine ioutmat4(a, l, m, n, nel, name6, ie)
Definition: hmholtz.f:1675
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
vlsc32
function vlsc32(r, b, m, n)
Definition: hmholtz.f:848
cggo_dg
subroutine cggo_dg(x, f, h1, h2, binv, mask, name, tin, maxit)
Definition: hmholtz.f:1494
set_fdm_prec_h1b
subroutine set_fdm_prec_h1b(d, h1, h2, nel)
Definition: hmholtz.f:1274
hmholtz_dg
subroutine hmholtz_dg(name, u, rhs, h1, h2, mask, tol, maxit)
Definition: hmholtz.f:1463
face_diff
subroutine face_diff(u, d, gsh_loc, w)
Definition: hmholtz.f:1772
set_fdm_prec_h1a_gen
subroutine set_fdm_prec_h1a_gen
Definition: hmholtz.f:1028
ioutfld
subroutine ioutfld(a, m, n, nel, name6, ie)
Definition: hmholtz.f:1695
set_fdm_prec_h1a
subroutine set_fdm_prec_h1a
Definition: hmholtz.f:1360
calc
subroutine calc(diag, upper, d, e, n, dmax, dmin)
Definition: hmholtz.f:858
setprec_dg
subroutine setprec_dg(d, h1, h2, imsh, isd)
Definition: hmholtz.f:1799
rescale_abhalf
subroutine rescale_abhalf(a, b, w, n)
Definition: hmholtz.f:1436
setprec
subroutine setprec(dpcm1, helm1, helm2, imsh, isd)
Definition: hmholtz.f:381
h1mg_solve
subroutine h1mg_solve(z, rhs, if_hybrid)
Definition: hsmg.f:1856
cadd
subroutine cadd(a, const, n)
Definition: math.f:326
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
glmin
function glmin(a, n)
Definition: math.f:973
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
glsc2
function glsc2(x, y, n)
Definition: math.f:794
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
vlamax
real function vlamax(vec, n)
Definition: math.f:406
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
addcol4
subroutine addcol4(a, b, c, d, n)
Definition: math.f:142
add3
subroutine add3(a, b, c, n)
Definition: math.f:644
col4
subroutine col4(a, b, c, d, n)
Definition: math.f:120
vlsum
real function vlsum(vec, n)
Definition: math.f:417
vsqrt
subroutine vsqrt(A, N)
Definition: math.f:41
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
glsum
function glsum(x, n)
Definition: math.f:861
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
glamax
real function glamax(a, n)
Definition: math.f:874
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
rone
subroutine rone(a, n)
Definition: math.f:230
invcol1
subroutine invcol1(a, n)
Definition: math.f:62
add2s1
subroutine add2s1(a, b, c1, n)
Definition: math.f:678
vlmin
real function vlmin(vec, n)
Definition: math.f:357
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
mxma
subroutine mxma(a, n1, b, n2, c, n3)
Definition: mxm_std.f:4118
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
ortho
subroutine ortho(respr)
Definition: navier1.f:224
vlsc3
function vlsc3(X, Y, B, N)
Definition: navier4.f:324
crs_solve_h1
subroutine crs_solve_h1(uf, vf)
Definition: navier8.f:1341
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342