KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier4.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2 c
3 c TO DO: Need to monitor dssum when using DG. (11/29/15, pff)
4 c
5 c Specifically - line 837 in navier4.f
6 c - line 537 in hmholtz.f
7 c
8 c-----------------------------------------------------------------------
9  subroutine setrhs(p,h1,h2,h2inv)
10 C
11 C Project rhs onto best fit in the "E" norm.
12 C
13  include 'SIZE'
14  include 'INPUT'
15  include 'MASS'
16  include 'SOLN'
17  include 'TSTEP'
18 C
19  REAL P (LX2,LY2,LZ2,LELV)
20  REAL H1 (LX1,LY1,LZ1,LELV)
21  REAL H2 (LX1,LY1,LZ1,LELV)
22  REAL H2INV(LX1,LY1,LZ1,LELV)
23 C
24  logical ifdump
25  save ifdump
26  data ifdump /.false./
27 C
28  parameter(ltot2=lx2*ly2*lz2*lelv)
29  COMMON /orthov/ rhs(ltot2,mxprev)
30  COMMON /orthox/ pbar(ltot2),pnew(ltot2)
31  COMMON /orthos/ alpha(mxprev), work(mxprev), alphan, dtlast
32  COMMON /orthoi/ nprev,mprev
33  REAL ALPHA,WORK
34 C
35 C
36  integer icalld
37  save icalld
38  data icalld/0/
39 C
40 C First call, we have no vectors to orthogonalize against.
41  IF (icalld.EQ.0) THEN
42  icalld=icalld+1
43  nprev=0
44  mprev=param(93)
45  mprev=min(mprev,mxprev)
46  endif
47 C
48 C Diag to see how much reduction in the residual is attained.
49 C
50  ntot2 = lx2*ly2*lz2*nelv
51  alpha1 = glsc3(p,p,bm2inv,ntot2)
52  if (alpha1.gt.0) alpha1 = sqrt(alpha1/volvm2)
53 C
54 C Update rhs's if E-matrix has changed
55 C
56  CALL updrhse(p,h1,h2,h2inv,ierr)
57  if (ierr.eq.1) nprev=0
58 C
59 C Perform Gram-Schmidt for previous rhs's.
60 C
61  DO 10 i=1,nprev
62  alpha(i) = vlsc2(p,rhs(1,i),ntot2)
63  10 CONTINUE
64 C
65  IF (nprev.GT.0) CALL gop(alpha,work,'+ ',nprev)
66 C
67  CALL rzero(pbar,ntot2)
68  DO 20 i=1,nprev
69  alphas = alpha(i)
70  CALL add2s2(pbar,rhs(1,i),alphas,ntot2)
71  20 CONTINUE
72 C
73  if (nprev.gt.0) then
74  intetype = 1
75  CALL cdabdtp(pnew,pbar,h1,h2,h2inv,intetype)
76  CALL sub2 (p,pnew,ntot2)
77 C ................................................................
78 C Diag.
79  alpha2 = glsc3(p,p,bm2inv,ntot2)
80  if (alpha2.gt.0) alpha2 = sqrt(alpha2/volvm2)
81  ratio = alpha1/alpha2
82  n10=min(10,nprev)
83 c IF (NIO.EQ.0) WRITE(6,11)ISTEP,Nprev,(ALPHA(I),I=1,N10)
84 c IF (NIO.EQ.0) WRITE(6,12) ISTEP,nprev,ALPHA1,ALPHA2,ratio
85  11 FORMAT(2i5,' alpha:',1p10e12.4)
86  12 FORMAT(i6,i4,1p3e12.4,' alph12')
87 C
88 c if (alpha1.gt.0.001 .and. .not.ifdump) then
89 c IF (NID.EQ.0) WRITE(6,*) 'alph1 large ... '
90 c if (istep.gt.10) then
91 c IF (NID.EQ.0) WRITE(6,*) ' ... dumping'
92 c call prepost (.true.,' ')
93 c ifdump = .true.
94 c else
95 c IF (NID.EQ.0) WRITE(6,*) ' ... doing nothing'
96 c endif
97 c endif
98 c if (alpha1.gt.0.1.and.istep.gt.10) then
99 c IF (NID.EQ.0) WRITE(6,*) 'alph1 too large ... aborting'
100 c call prepost (.true.,' ')
101 c call exitt
102 c endif
103 C ................................................................
104  endif
105 C
106 C
107  return
108  end
109 c-----------------------------------------------------------------------
110  subroutine gensoln(p,h1,h2,h2inv)
111 C
112 C Reconstruct the solution to the original problem by adding back
113 C the previous solutions
114 C know the soln.
115 C
116  include 'SIZE'
117  parameter(ltot2=lx2*ly2*lz2*lelv)
118  COMMON /orthov/ rhs(ltot2,mxprev)
119  COMMON /orthox/ pbar(ltot2),pnew(ltot2)
120  COMMON /orthos/ alpha(mxprev), work(mxprev), alphan, dtlast
121  COMMON /orthoi/ nprev,mprev
122  REAL ALPHA,WORK
123 
124  REAL P (LX2,LY2,LZ2,LELV)
125  REAL H1 (LX1,LY1,LZ1,LELV)
126  REAL H2 (LX1,LY1,LZ1,LELV)
127  REAL H2INV(LX1,LY1,LZ1,LELV)
128 C
129  ntot2=lx2*ly2*lz2*nelv
130 C
131 C First, save current solution
132 C
133  CALL copy (pnew,p,ntot2)
134 C
135 C Reconstruct solution
136 C
137  CALL add2(p,pbar,ntot2)
138 C
139 C Update the set of <p,rhs>
140 C
141  CALL updtset(p,h1,h2,h2inv,ierr)
142  if (ierr.eq.1) nprev = 0
143 c
144  return
145  end
146 c-----------------------------------------------------------------------
147  subroutine updtset(p,h1,h2,h2inv,IERR)
148 C
149 C Update the set of rhs's and the corresponding p-set:
150 C
151 C . Standard case is to add P_new, and RHS_new = E*P_new
152 C
153 C . However, when Nprev=Mprev (max. allowed), we throw out
154 C the old set, and set P_1 = P, RHS_1=E*P_1
155 C
156 C . Other schemes are possible, e.g., let's save a bunch of
157 C old vectors, perhaps chosen wisely via P.O.D.
158 C
159 C
160  include 'SIZE'
161  include 'INPUT'
162  include 'MASS'
163  parameter(ltot2=lx2*ly2*lz2*lelv)
164  COMMON /orthov/ rhs(ltot2,mxprev)
165  COMMON /orthox/ pbar(ltot2),pnew(ltot2)
166  COMMON /orthos/ alpha(mxprev), work(mxprev), alphan, dtlast
167  COMMON /orthoi/ nprev,mprev
168 
169  REAL ALPHA,WORK
170 
171  REAL P (LX2,LY2,LZ2,LELV)
172  REAL H1 (LX1,LY1,LZ1,LELV)
173  REAL H2 (LX1,LY1,LZ1,LELV)
174  REAL H2INV(LX1,LY1,LZ1,LELV)
175 C
176  ntot2=lx2*ly2*lz2*nelv
177 C
178  IF (nprev.EQ.mprev) THEN
179  CALL copy(pnew,p,ntot2)
180  nprev=0
181  endif
182 C
183 C Increment solution set
184  nprev = nprev+1
185 C
186  CALL copy (rhs(1,nprev),pnew,ntot2)
187 
188 C Orthogonalize rhs against previous rhs and normalize
189 
190  write(6,*) istep,nprev,' call econj2'
191  call econj (nprev,h1,h2,h2inv,ierr)
192 C call echeck(nprev,h1,h2,h2inv,intetype)
193 C
194 c Save last sol'n
195  CALL copy(pnew,p,ntot2)
196 C
197  return
198  end
199 c-----------------------------------------------------------------------
200  subroutine econj(kprev,h1,h2,h2inv,ierr)
201 C
202 C Orthogonalize the rhs wrt previous rhs's for which we already
203 C know the soln.
204 C
205  include 'SIZE'
206  include 'INPUT'
207  include 'MASS'
208  include 'SOLN'
209  include 'TSTEP'
210 C
211  REAL H1 (LX1,LY1,LZ1,LELV)
212  REAL H2 (LX1,LY1,LZ1,LELV)
213  REAL H2INV(LX1,LY1,LZ1,LELV)
214 C
215  parameter(ltot2=lx2*ly2*lz2*lelv)
216  COMMON /orthov/ rhs(ltot2,mxprev)
217  COMMON /orthox/ pbar(ltot2),pnew(ltot2),pbrr(ltot2)
218  COMMON /orthos/ alpha(mxprev), work(mxprev), alphan, dtlast
219  COMMON /orthoi/ nprev,mprev
220  REAL ALPHA,WORK
221  real ALPHAd
222 C
223 C
224  ierr = 0
225  ntot2 = lx2*ly2*lz2*nelv
226  intetype=1
227 C
228 C Gram Schmidt, w re-orthogonalization
229 C
230  npass=1
231  if (abs(param(105)).eq.2) npass=2
232  do ipass=1,npass
233 c
234  CALL cdabdtp(pbrr,rhs(1,kprev),h1,h2,h2inv,intetype)
235 C
236 C Compute part of the norm
237  alphad = glsc2(rhs(1,kprev),pbrr,ntot2)
238 C
239 C Gram-Schmidt
240  kprev1=kprev-1
241  DO 10 i=1,kprev1
242  alpha(i) = vlsc2(pbrr,rhs(1,i),ntot2)
243  10 CONTINUE
244  IF (kprev1.GT.0) CALL gop(alpha,work,'+ ',kprev1)
245 C
246  DO 20 i=1,kprev1
247  alpham = -alpha(i)
248  CALL add2s2(rhs(1,kprev),rhs(1,i),alpham,ntot2)
249  alphad = alphad - alpha(i)**2
250  20 CONTINUE
251  enddo
252 C
253 C .Normalize new element in P~
254 C
255  if (alphad.le.0.0) then
256  if (nio.eq.0)
257  $ write(6,*) .le.'ERROR: alphad 0 in econj',alphad,kprev
258  ierr = 1
259  return
260  endif
261  alphad = 1.0/sqrt(alphad)
262  alphan = alphad
263  CALL cmult(rhs(1,kprev),alphan,ntot2)
264 C
265  return
266  end
267 c-----------------------------------------------------------------------
268  subroutine chkptol
269 C--------------------------------------------------------------------
270 C
271 C Check pressure tolerance for transient case.
272 C
273 C pff 6/20/92
274 C This routine has been modified for diagnostic purposes only.
275 C It can be replaced with the standard nekton version.
276 C
277 C--------------------------------------------------------------------
278  include 'SIZE'
279  include 'SOLN'
280  include 'MASS'
281  include 'INPUT'
282  include 'TSTEP'
283  COMMON /ctolpr/ divex
284  COMMON /cprint/ ifprint
285  LOGICAL IFPRINT
286 C
287  COMMON /scruz/ divv(lx2,ly2,lz2,lelv)
288  $ , bdivv(lx2,ly2,lz2,lelv)
289 C
290  if (ifsplit) return
291  IF (param(102).eq.0.and.(tolpdf.NE.0. .OR. istep.LE.5)) return
292  five = 5.0
293  if (param(102).ne.0.0) five=param(102)
294 C
295  ntot2 = lx2*ly2*lz2*nelv
296  if (ifield.eq.1) then ! avo: sub arguments?
297  CALL opdiv (bdivv,vx,vy,vz)
298  else
299  CALL opdiv (bdivv,bx,by,bz)
300  endif
301  CALL col3 (divv,bdivv,bm2inv,ntot2)
302  dnorm = sqrt(glsc2(divv,bdivv,ntot2)/volvm2)
303 C
304  if (nio.eq.0) WRITE (6,*) istep,' DNORM, DIVEX',dnorm,divex
305 C
306 c IF (istep.gt.10.and.DNORM.GT.(1.01*DIVEX).AND.DIVEX.GT.0.) then
307 c if (DNORM.gt.1e-8) then
308 c if (nid.eq.0) WRITE(6,*) 'DNORM-DIVEX div. ... aborting'
309 c call prepost (.true.,' ')
310 c call exitt
311 c else
312 c if (nid.eq.0) WRITE(6,*) 'DNORM-DIVEX div. ... small'
313 c endif
314 c endif
315 C
316 c IF (DNORM.GT.(1.2*DIVEX).AND.DIVEX.GT.0.) TOLPDF = 5.*DNORM
317  IF (istep.gt.5.and.tolpdf.eq.0.0.and.
318  $ dnorm.GT.(1.2*divex).AND.divex.GT.0.)
319  $ tolpdf = five*dnorm
320 C
321  return
322  end
323  FUNCTION vlsc3(X,Y,B,N)
324 C
325 C local inner product, with weight
326 C
327  dimension x(1),y(1),b(1)
328  REAL dt
329 C
330  include 'OPCTR'
331 C
332 #ifdef TIMER
333 C
334  if (isclld.eq.0) then
335  isclld=1
336  nrout=nrout+1
337  myrout=nrout
338  rname(myrout) = 'VLSC3 '
339  endif
340  isbcnt = 3*n
341  dct(myrout) = dct(myrout) + dfloat(isbcnt)
342  ncall(myrout) = ncall(myrout) + 1
343  dcount = dcount + dfloat(isbcnt)
344 #endif
345 C
346  dt = 0.0
347  DO 10 i=1,n
348  t = x(i)*y(i)*b(i)
349  dt = dt+t
350  10 CONTINUE
351  t=dt
352  vlsc3 = t
353  return
354  end
355 c-----------------------------------------------------------------------
356  subroutine updrhse(p,h1,h2,h2inv,ierr)
357 C
358 C Update rhs's if E-matrix has changed
359 C
360 C
361  include 'SIZE'
362  include 'INPUT'
363  include 'MASS'
364  include 'TSTEP'
365 C
366  parameter(ltot2=lx2*ly2*lz2*lelv)
367  COMMON /orthov/ rhs(ltot2,mxprev)
368  COMMON /orthox/ pbar(ltot2),pnew(ltot2)
369  COMMON /orthos/ alpha(mxprev), work(mxprev), alphan, dtlast
370  COMMON /orthoi/ nprev,mprev
371  COMMON /orthol/ ifnewe
372  REAL ALPHA,WORK
373  LOGICAL IFNEWE
374 C
375 C
376  REAL P (LX2,LY2,LZ2,LELV)
377  REAL H1 (LX1,LY1,LZ1,LELV)
378  REAL H2 (LX1,LY1,LZ1,LELV)
379  REAL H2INV(LX1,LY1,LZ1,LELV)
380 C
381  integer icalld
382  save icalld
383  data icalld/0/
384 
385  ntot2=lx2*ly2*lz2*nelv
386 
387 
388 C First, we have to decide if the E matrix has changed.
389 
390 
391  if (icalld.eq.0) then
392  icalld=1
393  dtlast=dt
394  endif
395 
396  ifnewe=.false.
397  if (ifmvbd) then
398  ifnewe=.true.
399  call invers2(bm2inv,bm2,ntot2)
400  endif
401  if (dtlast.ne.dt) then
402  ifnewe=.true.
403  dtlast=dt
404  endif
405 c if (ifnewe.and.nio.eq.0) write(6,*) istep,'reorthogo:',nprev
406 
407 
408 C
409 C Next, we reconstruct a new rhs set.
410 C
411  if (ifnewe) then
412 c
413 c new idea...
414 c if (nprev.gt.0) nprev=1
415 c call copy(rhs,pnew,ntot2)
416 c
417  nprevt = nprev
418  DO 100 iprev=1,nprevt
419 C Orthogonalize this rhs w.r.t. previous rhs's
420  call econj (iprev,h1,h2,h2inv,ierr)
421  if (ierr.eq.1) then
422  if (nio.eq.0) write(6,*) istep,ierr,' econj error'
423  nprev = 0
424  return
425  endif
426  100 CONTINUE
427 C
428  endif
429 C
430  return
431  end
432 c-----------------------------------------------------------------------
433  subroutine hconj(approx,k,h1,h2,vml,vmk,ws,name4,ierr)
434 c
435 c Orthonormalize the kth vector against vector set
436 c
437  include 'SIZE'
438  include 'INPUT'
439  include 'MASS'
440  include 'SOLN'
441  include 'PARALLEL'
442  include 'TSTEP'
443 c
444  parameter(lt=lx1*ly1*lz1*lelt)
445  real approx(lt,0:1),h1(1),h2(1),vml(1),vmk(1),ws(1)
446  character*4 name4
447 c
448  ierr=0
449  ntot=lx1*ly1*lz1*nelfld(ifield)
450 c
451  call axhelm (approx(1,0),approx(1,k),h1,h2,1,1)
452  call col2 (approx(1,0),vmk,ntot)
453 c
454 c Compute part of the norm (Note: a(0) already scaled by vml)
455 c
456  alpha = glsc2(approx(1,0),approx(1,k),ntot)
457  alph1 = alpha
458 c
459 c Gram-Schmidt
460 c
461  km1=k-1
462  do i=1,km1
463  ws(i) = vlsc2(approx(1,0),approx(1,i),ntot)
464  enddo
465  if (km1.gt.0) call gop(ws,ws(k),'+ ',km1)
466 c
467  do i=1,km1
468  alpham = -ws(i)
469  call add2s2(approx(1,k),approx(1,i),alpham,ntot)
470  alpha = alpha - ws(i)**2
471  enddo
472 c
473 c .Normalize new element in approximation space
474 c
475  eps = 1.e-7
476  if (wdsize.eq.8) eps = 1.e-15
477  ratio = alpha/alph1
478 c
479  if (ratio.le.0) then
480  ierr=1
481  if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
482  12 format(i6,1x,a4,' alpha b4 sqrt:',i4,1p2e12.4)
483  elseif (ratio.le.eps) then
484  ierr=2
485  if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
486  else
487  ierr=0
488  alpha = 1.0/sqrt(alpha)
489  call cmult(approx(1,k),alpha,ntot)
490  endif
491 c
492  if (ierr.ne.0) then
493  call axhelm (approx(1,0),approx(1,k),h1,h2,1,1)
494  call col2 (approx(1,0),vmk,ntot)
495 c
496 c Compute part of the norm (Note: a(0) already scaled by vml)
497 c
498  alpha = glsc2(approx(1,0),approx(1,k),ntot)
499  if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
500  if (alpha.le.0) then
501  ierr=3
502  if (nio.eq.0) write(6,12) istep,name4,k,alpha,alph1
503  return
504  endif
505  alpha = 1.0/sqrt(alpha)
506  call cmult(approx(1,k),alpha,ntot)
507  ierr = 0
508  endif
509 c
510  return
511  end
512 c-----------------------------------------------------------------------
513  subroutine hmhzpf(name,u,r,h1,h2,mask,mult,imesh,tli,maxit,isd,bi)
514  include 'SIZE'
515  include 'INPUT'
516  include 'MASS'
517  include 'FDMH1'
518  include 'CTIMER'
519 c
520  CHARACTER*4 NAME
521  REAL U (LX1,LY1,LZ1,1)
522  REAL R (LX1,LY1,LZ1,1)
523  REAL H1 (LX1,LY1,LZ1,1)
524  REAL H2 (LX1,LY1,LZ1,1)
525  REAL MASK (LX1,LY1,LZ1,1)
526  REAL MULT (LX1,LY1,LZ1,1)
527  REAL bi (LX1,LY1,LZ1,1)
528  COMMON /ctmp0/ w1(lx1,ly1,lz1,lelt)
529  $ , w2(lx1,ly1,lz1,lelt)
530 c
531  etime1=dnekclock()
532 c
533  IF (imesh.EQ.1) ntot = lx1*ly1*lz1*nelv
534  IF (imesh.EQ.2) ntot = lx1*ly1*lz1*nelt
535 c
536  tol = tli
537  if (param(22).ne.0) tol = abs(param(22))
538  CALL chktcg1 (tol,r,h1,h2,mask,mult,imesh,isd)
539 c
540 c
541 c Set flags for overlapping Schwarz preconditioner (pff 11/12/98)
542 c
543  kfldfdm = -1
544 c if (name.eq.'TEMP') kfldfdm = 0
545 c if (name.eq.'VELX') kfldfdm = 1
546 c if (name.eq.'VELY') kfldfdm = 2
547 c if (name.eq.'VELZ') kfldfdm = 3
548  if (name.eq.'PRES') kfldfdm = ldim+1
549 
550  if (ifdg) then
551  call cggo_dg (u,r,h1,h2,bi,mask,name,tol,maxit)
552  else
553  call cggo
554  $ (u,r,h1,h2,mask,mult,imesh,tol,maxit,isd,bi,name)
555  endif
556  thmhz=thmhz+(dnekclock()-etime1)
557 c
558 c
559  return
560  end
561 c-----------------------------------------------------------------------
562  subroutine hsolve(name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd
563  $ ,approx,napprox,bi)
564 c
565 c Either std. Helmholtz solve, or a projection + Helmholtz solve
566 c
567  include 'SIZE'
568  include 'TOTAL'
569  include 'CTIMER'
570 c
571  CHARACTER*4 NAME
572  REAL U (LX1,LY1,LZ1,1)
573  REAL R (LX1,LY1,LZ1,1)
574  REAL H1 (LX1,LY1,LZ1,1)
575  REAL H2 (LX1,LY1,LZ1,1)
576  REAL vmk (LX1,LY1,LZ1,1)
577  REAL vml (LX1,LY1,LZ1,1)
578  REAL bi (LX1,LY1,LZ1,1)
579  REAL approx (1)
580  integer napprox(1)
581  common /ctmp2/ w1(lx1,ly1,lz1,lelt)
582  common /ctmp3/ w2(2+2*mxprev)
583 
584  logical ifstdh
585  character*4 cname
586  character*6 name6
587 
588  logical ifwt,ifvec
589 
590  call chcopy(cname,name,4)
591  call capit (cname,4)
592 
593  ifstdh = .true. ! default is no projection
594  if (ifprojfld(ifield)) then
595  ifstdh = .false.
596  endif
597 
598  p945 = param(94)
599  if (cname.eq.'PRES') then
600  ifstdh = .false.
601  p945 = param(95)
602  endif
603 
604  if (ifield.gt.ldimt_proj+1) ifstdh = .true.
605  if (param(93).eq.0) ifstdh = .true.
606  if (p945.eq.0) ifstdh = .true.
607  if (istep.lt.p945) ifstdh = .true.
608 
609  if (ifstdh) then
610  call hmholtz(name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd)
611  else
612 
613  n = lx1*ly1*lz1*nelfld(ifield)
614 
615  call col2 (r,vmk,n)
616  call dssum (r,lx1,ly1,lz1)
617 
618  call blank (name6,6)
619  call chcopy(name6,name,4)
620  ifwt = .true.
621  ifvec = .false.
622 
623  call project1
624  $ (r,n,approx,napprox,h1,h2,vmk,vml,ifwt,ifvec,name6)
625 
626  call hmhzpf (name,u,r,h1,h2,vmk,vml,imsh,tol,maxit,isd,bi)
627 
628  call project2
629  $ (u,n,approx,napprox,h1,h2,vmk,vml,ifwt,ifvec,name6)
630 
631  endif
632 
633  return
634  end
635 c-----------------------------------------------------------------------
636  subroutine project1(b,n,rvar,ivar,h1,h2,msk,w,ifwt,ifvec,name6)
637 
638 c 1. Compute the projection of x onto X
639 
640 c 2. Re-orthogonalize the X basis set and corresponding B=A*X
641 c vectors if A has changed.
642 
643 c Output: b = b - projection of b onto B
644 
645 c Input: n = length of field (or multifields, when ifvec=true)
646 c rvar = real array of field values, including old h1,h2, etc.
647 c ivar = integer array of pointers, etc.
648 c h1 = current h1, for Axhelm(.,.,h1,h2,...)
649 c h2 = current h2
650 c msk = mask for Dirichlet BCs
651 c w = weight for inner products (typ. w=vmult, tmult, etc.)
652 c ifwt = use weighted inner products when ifwt=.true.
653 c ifvec = are x and b vectors, or scalar fields?
654 c name6 = discriminator for action of A*x
655 
656 c The idea here is to have one pair of projection routines for
657 c constructing the new rhs (project1) and reconstructing the new
658 c solution (x = xbar + dx) plus updating the approximation space.
659 c The latter functions are done in project2.
660 c
661 c The approximation space X and corresponding right-hand sides,
662 c B := A*X are stored in rvar, as well as h1old and h2old and a
663 c couple of other auxiliary arrays.
664 
665 c In this new code, we retain both X and B=A*X and we re-orthogonalize
666 c at each timestep (with no extra matrix-vector products, but O(nm)
667 c work. The idea is to retain fresh vectors by injecting the most
668 c recent solution and pushing the oldest off the stack, hopefully
669 c keeping the number of vectors, m, small.
670 
671 
672  include 'SIZE' ! For nid/nio
673  include 'TSTEP' ! For istep
674  include 'CTIMER'
675 
676  real b(n),rvar(n,1),h1(n),h2(n),w(n),msk(n)
677  integer ivar(1)
678  character*6 name6
679  logical ifwt,ifvec
680 
681  etime0 = dnekclock()
682 
683  nn = n
684  if (ifvec) nn = n*ldim
685 
686  call proj_get_ivar
687  $ (m,mmx,ixb,ibb,ix,ib,ih1,ih2,ivar,n,ifvec,name6)
688 
689  if (m.le.0) return
690 
691  ireset=iproj_chk(rvar(ih1,1),rvar(ih2,1),h1,h2,n) ! Updated matrix?
692 
693  bb4 = glsc3(b,w,b,n)
694  bb4 = sqrt(bb4)
695 
696 
697 c Re-orthogonalize basis set w.r.t. new vectors if space has changed.
698 
699  if (ireset.eq.1) then
700 
701  do j=0,m-1 ! First, set B := A*X
702  jb = ib+j*nn !Iterate through xs and bs
703  jx = ix+j*nn
704  call proj_matvec (rvar(jb,1),rvar(jx,1),n,h1,h2,msk,name6)
705  enddo
706 
707  if (nio.eq.0 .and. loglevel.gt.2)
708  $ write(6,'(13x,A)') 'Reorthogonalize Basis'
709 
710  call proj_ortho_full
711  $ (rvar(ix,1),rvar(ib,1),n,m,w,ifwt,ifvec,name6)
712 
713  ivar(2) = m ! Update number of saved vectors
714 
715  endif
716 
717 c ixb is pointer to xbar, ibb is pointer to bbar := A*xbar
718 
719  call project1_a(rvar(ixb,1),rvar(ibb,1),b,rvar(ix,1),rvar(ib,1)
720  $ ,n,m,w,ifwt,ifvec)
721 
722  baf = glsc3(b,w,b,n)
723  baf = sqrt(baf)
724  ratio = bb4/baf
725 
726  tproj = tproj + dnekclock() - etime0
727 
728  if (nio.eq.0) write(6,1) istep,' Project ' // name6,
729  & baf,bb4,ratio,m,mmx
730  1 format(i11,a,6x,1p3e13.4,i4,i4)
731 
732  return
733  end
734 c-----------------------------------------------------------------------
735  subroutine project1_a(xbar,bbar,b,xx,bb,n,m,w,ifwt,ifvec)
736 
737 c xbar is best fit in xx, bbar = A*xbar
738 c b <-- b - bbar
739 
740  include 'SIZE'
741  include 'TSTEP'
742  include 'INPUT'
743  include 'PARALLEL'
744  parameter(lt=lx1*ly1*lz1*lelv)
745  real xbar(n),bbar(n),b(n),xx(n,m),bb(n,m),w(n)
746  logical ifwt,ifvec
747  integer e,eg
748  parameter(lxyz=lx1*ly1*lz1)
749  real tt(lxyz,lelt)
750 
751  real work(mxprev),alpha(mxprev)
752 
753  if (m.le.0) return
754 
755  !First round of CGS
756  do k = 1, m
757  if(ifwt) then
758  alpha(k) = vlsc3(xx(1,k),w,b,n)
759  else
760  alpha(k) = vlsc2(xx(1,k),b,n)
761  endif
762  enddo
763  !First one outside loop to avoid zeroing xbar and bbar
764  call gop(alpha,work,'+ ',m)
765  call cmult2(xbar,xx(1,1),alpha(1),n)
766  call cmult2(bbar,bb(1,1),alpha(1),n)
767  call add2s2(b,bb(1,1),-alpha(1),n)
768  do k = 2,m
769  call add2s2(xbar,xx(1,k),alpha(k),n)
770  call add2s2(bbar,bb(1,k),alpha(k),n)
771  call add2s2(b,bb(1,k),-alpha(k),n)
772  enddo
773  !Second round of CGS
774  do k = 1, m
775  if(ifwt) then
776  alpha(k) = vlsc3(xx(1,k),w,b,n)
777  else
778  alpha(k) = vlsc2(xx(1,k),b,n)
779  endif
780  enddo
781  call gop(alpha,work,'+ ',m)
782  do k = 1,m
783  call add2s2(xbar,xx(1,k),alpha(k),n)
784  call add2s2(bbar,bb(1,k),alpha(k),n)
785  call add2s2(b,bb(1,k),-alpha(k),n)
786  enddo
787 
788  return
789  end
790 c-----------------------------------------------------------------------
791  function iproj_chk(h1old,h2old,h1,h2,n)
792  include 'SIZE'
793  include 'TOTAL'
794 
795 c Matrix has changed if h1/h2 differ from old values
796 
797  real h1(n),h2(n),h1old(n),h2old(n)
798 
799  iproj_chk = 0
800 
801  if (ifmvbd) then
802  iproj_chk = 1
803  return
804  endif
805 
806  dh1 = 0.
807  dh2 = 0.
808  do i=1,n
809  dh1 = max(dh1,abs(h1(i)-h1old(i)))
810  dh2 = max(dh2,abs(h2(i)-h2old(i)))
811  enddo
812  dh = max(dh1,dh2)
813  dh = glmax(dh,1) ! Max across all processors
814 
815  if (dh.gt.0) then
816 
817  call copy(h1old,h1,n) ! Save old h1 / h2 values
818  call copy(h2old,h2,n)
819 
820  iproj_chk = 1 ! Force re-orthogonalization of basis
821 
822  endif
823 
824  return
825  end
826 c-----------------------------------------------------------------------
827  subroutine proj_matvec(b,x,n,h1,h2,msk,name6)
828  include 'SIZE'
829  include 'TOTAL'
830  real b(n),x(n),h1(n),h2(n),msk(n)
831  character*6 name6
832 
833 c This is the default matvec for nekcem.
834 
835 c The code can later be updated to support different matvec
836 c implementations, which would be discriminated by the character
837 c string "name6"
838 
839  isd = 1 ! This probably won't work for axisymmetric
840  imsh = 1
841  if (iftmsh(ifield)) imsh=2
842 
843  call axhelm (b,x,h1,h2,imsh,isd) ! b = A x
844  call dssum (b,lx1,ly1,lz1)
845  call col2 (b,msk,n)
846 
847  return
848  end
849 c-----------------------------------------------------------------------
850  !O(nm) method for updating projection space
851  !See James Lotte's note or Nicholas Christensen's master's thesis
852  subroutine proj_ortho(xx,bb,n,m,w,ifwt,ifvec,name6)
853 
854  include 'SIZE' ! nio
855  include 'TSTEP' ! istep
856  include 'PARALLEL' ! wdsize
857 
858  real xx(n,1), bb(n,1), w(n)
859  character*6 name6
860  logical ifwt, ifvec
861  real tol, nrm, scl1, scl2, c, s
862  real work(mxprev), alpha(mxprev), beta(mxprev)
863  integer h
864 
865  if(m.le.0) return !No vectors to ortho-normalize
866 
867  ! AX = B
868  ! Calculate dx, db: dx = x-XX^Tb, db=b-BX^Tb
869 
870  do k = 1, m !First round CGS
871  if(ifwt) then
872  alpha(k) = 0.5*(vlsc3(xx(1,k),w,bb(1,m),n)
873  $ + vlsc3(bb(1,k),w,xx(1,m),n))
874  else
875  alpha(k) = 0.5*(vlsc2(xx(1,k),bb(1,m),n)
876  $ + vlsc2(bb(1,k),xx(1,m),n))
877  endif
878  enddo
879  call gop(alpha,work,'+ ',m)
880  nrm = sqrt(alpha(m)) !Calculate A-norm of new vector
881  do k = 1,m-1
882  call add2s2(xx(1,m),xx(1,k),-alpha(k),n)
883  call add2s2(bb(1,m),bb(1,k),-alpha(k),n)
884  enddo
885 
886  do k = 1, m-1 !Second round CGS
887  if(ifwt) then
888  beta(k) = 0.5*(vlsc3(xx(1,k),w,bb(1,m),n)
889  $ + vlsc3(bb(1,k),w,xx(1,m),n))
890  else
891  beta(k) = 0.5*(vlsc2(xx(1,k),bb(1,m),n)
892  $ + vlsc2(bb(1,k),xx(1,m),n))
893  endif
894  enddo
895  call gop(beta,work,'+ ',m-1)
896  do k = 1,m-1
897  call add2s2(xx(1,m),xx(1,k),-beta(k),n)
898  call add2s2(bb(1,m),bb(1,k),-beta(k),n)
899  !While we're at it,
900  !Sum weights from each round to get the total alpha
901  alpha(k) = alpha(k) + beta(k)
902  enddo
903 
904  !Calculate A-norm of newest solution
905  if(ifwt) then
906  alpha(m) = glsc3(xx(1,m), w, bb(1,m), n)
907  else
908  alpha(m) = glsc2(xx(1,m), bb(1,m), n)
909  endif
910  alpha(m) = sqrt(alpha(m))
911  !dx and db now stored in last column of xx and bb
912 
913 c Set tolerance for linear independence
914  tol = 1.e-7
915  if (wdsize.eq.4) tol=1.e-3
916 
917 c Check for linear independence.
918  if(alpha(m).gt.tol*nrm) then !New vector is linearly independent
919 
920  !Normalize dx and db
921  scl1 = 1.0/alpha(m)
922  call cmult(xx(1,m), scl1, n)
923  call cmult(bb(1,m), scl1, n)
924 
925  !We want to throw away the oldest information
926  !The below propagates newest information to first vector.
927  !This will make the first vector a scalar
928  !multiple of x.
929  do k = m, 2, -1
930  h = k - 1
931  call givens_rotation(alpha(h),alpha(k),c,s,nrm)
932  alpha(h) = nrm
933  do i = 1, n !Apply rotation to xx and bb
934  scl1 = c*xx(i,h) + s*xx(i,k)
935  xx(i,k) = -s*xx(i,h) + c*xx(i,k)
936  xx(i,h) = scl1
937  scl2 = c*bb(i,h) + s*bb(i,k)
938  bb(i,k) = -s*bb(i,h) + c*bb(i,k)
939  bb(i,h) = scl2
940  enddo
941  enddo
942 
943  else !New vector is not linearly independent, forget about it
944  k = m !location of rank deficient column
945 
946  if (nio.eq.0 .and. loglevel.gt.2)
947  $ write(6,1) istep,k,m,name6,alpha(m),tol
948 
949  1 format(i9,'proj_ortho: ',2i4,1x,a6,
950  $ ' Detect rank deficiency:',1p2e12.4)
951 
952  m = m - 1 !Remove column
953  endif
954 
955  return
956  end
957 c-----------------------------------------------------------------------
958  !Function to switch between mgs and cgs2 full reorthogonalization
959  subroutine proj_ortho_full(xx,bb,n,m,w,ifwt,ifvec,name6)
960 
961  include 'SIZE'
962 
963  real xx(n,1),bb(n,1),w(n)
964  character*6 name6
965  logical ifwt,ifvec
966  integer flag(mxprev)
967 
968  !call proj_ortho_full_mgs(xx,bb,n,m,w,ifwt,ifvec,name6)
969  call proj_ortho_full_cgs2(xx,bb,n,m,w,ifwt,ifvec,name6)
970 
971  return
972  end
973 c-----------------------------------------------------------------------
974 c Full MGS reorthogonalization
975  subroutine proj_ortho_full_mgs(xx,bb,n,m,w,ifwt,ifvec,name6)
976 
977  include 'SIZE' ! nio
978  include 'TSTEP' ! istep
979  include 'PARALLEL' ! wdsize
980 
981  real xx(n,1),bb(n,1),w(n)
982  character*6 name6
983  logical ifwt,ifvec
984  integer flag(mxprev)
985  real normk,normp,alpha
986 
987  if (m.le.0) return
988 
989  if ( ifwt) alpha = glsc3(xx(1,m),w,bb(1,m),n)
990  if (.not. ifwt) alpha = glsc2(xx(1,m),bb(1,m),n)
991  if (alpha.eq.0) return
992 
993  scale = 1./sqrt(alpha)
994  call cmult(xx(1,m),scale,n)
995  call cmult(bb(1,m),scale,n)
996  flag(m) = 1
997 
998  do k=m-1,1,-1 ! Reorthogonalize, starting with latest solution
999 
1000  if ( ifwt) normk = glsc3(xx(1,k),w,bb(1,k),n)
1001  if (.not. ifwt) normk = glsc2(xx(1,k),bb(1,k),n)
1002  normk=sqrt(normk)
1003 
1004  do j=m,k+1,-1 ! Modified GS
1005  if(flag(j).eq.1) then
1006  alpha = 0.
1007  if (ifwt) then
1008  alpha = alpha + .5*(vlsc3(xx(1,j),w,bb(1,k),n)
1009  $ + vlsc3(bb(1,j),w,xx(1,k),n))
1010  else
1011  alpha = alpha + .5*(vlsc2(xx(1,j),bb(1,k),n)
1012  $ + vlsc2(bb(1,j),xx(1,k),n))
1013  endif
1014  scale = -glsum(alpha,1)
1015  call add2s2(xx(1,k),xx(1,j),scale,n)
1016  call add2s2(bb(1,k),bb(1,j),scale,n)
1017  endif
1018  enddo
1019  if ( ifwt) normp = glsc3(xx(1,k),w,bb(1,k),n)
1020  if (.not. ifwt) normp = glsc2(xx(1,k),bb(1,k),n)
1021  if (normp.gt.0.0) normp=sqrt(normp)
1022 
1023  tol = 1.e-7
1024  if (wdsize.eq.4) tol=1.e-3
1025 
1026  if (normp.gt.tol*normk) then ! linearly independent vectors
1027  scale = 1./normp
1028  call cmult(xx(1,k),scale,n)
1029  call cmult(bb(1,k),scale,n)
1030  flag(k) = 1
1031 c if (nio.eq.0) write(6,2) istep,k,m,name6,normp,normk
1032 c 2 format(i9,'proj_ortho: ',2i4,1x,a6,' project ok.'
1033 c $ ,1p2e12.4)
1034  else
1035  flag(k) = 0
1036 c if (nio.eq.0) write(6,1) istep,k,m,name6,normp,normk
1037 c 1 format(i9,'proj_ortho: ',2i4,1x,a6,' Detect rank deficiency:'
1038 c $ ,1p2e12.4)
1039  endif
1040 
1041  enddo
1042 
1043  k=0
1044  do j=1,m
1045  if (flag(j).eq.1) then
1046  k=k+1
1047  if (k.lt.j) then
1048  call copy(xx(1,k),xx(1,j),n)
1049  call copy(bb(1,k),bb(1,j),n)
1050  endif
1051  endif
1052  enddo
1053  m = k
1054 
1055  return
1056  end
1057 c-----------------------------------------------------------------------
1058  !CGS2 version of full reorthogonalization, possibly more stable in
1059  !certain instances. Much faster for large m.
1060  subroutine proj_ortho_full_cgs2(xx,bb,n,m,w,ifwt,ifvec,name6)
1061 
1062  include 'SIZE' ! nio
1063  include 'TSTEP' ! istep
1064  include 'PARALLEL' ! wdsize
1065 
1066  real xx(n,1),bb(n,1),w(n)
1067  character*6 name6
1068  logical ifwt,ifvec
1069  integer flag(mxprev)
1070  real normk,normp,alpha(mxprev),work(mxprev),scl1,tol
1071 
1072  if (m.le.0) return
1073 
1074  tol = 1.e-7
1075  if (wdsize.eq.4) tol=1.e-3
1076 
1077  do i = 1, 2 !Do this twice for CGS2
1078 
1079  do k = m, 1, -1
1080  do j = m, k, -1
1081  alpha(j) = 0.0
1082  if(ifwt) then
1083  alpha(j) = .5*(vlsc3(xx(1,j),w,bb(1,k),n)
1084  $ + vlsc3(bb(1,j),w,xx(1,k),n))
1085  else
1086  alpha(j) = .5*(vlsc2(xx(1,j),bb(1,k),n)
1087  $ + vlsc2(bb(1,j),xx(1,k),n))
1088  endif
1089  enddo
1090  call gop(alpha(k), work,'+ ',(m - k) + 1)
1091  do j = m, k+1, -1
1092  call add2s2(xx(1,k),xx(1,j),-alpha(j),n)
1093  call add2s2(bb(1,k),bb(1,j),-alpha(j),n)
1094  enddo
1095  normp = sqrt(alpha(k))
1096  if(ifwt) then
1097  normk = glsc3(xx(1,k),w,bb(1,k),n)
1098  else
1099  normk = glsc2(xx(1,k),bb(1,k),n)
1100  endif
1101  normk = sqrt(normk)
1102  if(normk.gt.tol*normp) then
1103  scl1 = 1.0/normk
1104  call cmult(xx(1,k), scl1, n)
1105  call cmult(bb(1,k), scl1, n)
1106  flag(k) = 1
1107  else
1108  flag(k) = 0
1109  endif
1110  enddo
1111 
1112  enddo
1113 
1114  k=0
1115  do j=1,m
1116  if (flag(j).eq.1) then
1117  k=k+1
1118  if (k.lt.j) then
1119  call copy(xx(1,k),xx(1,j),n)
1120  call copy(bb(1,k),bb(1,j),n)
1121  endif
1122  endif
1123  enddo
1124  m = k
1125 
1126  return
1127  end
1128 c-----------------------------------------------------------------------
1129  subroutine project2(x,n,rvar,ivar,h1,h2,msk,w,ifwt,ifvec,name6)
1130 
1131  include 'CTIMER'
1132 
1133  real x(n),b(n),rvar(n,1),h1(n),h2(n),w(n),msk(n)
1134  integer ivar(1)
1135  character*6 name6
1136  logical ifwt,ifvec
1137 
1138  etime0 = dnekclock()
1139 
1140  call proj_get_ivar(m,mmx,ixb,ibb,ix,ib,ih1,ih2,
1141  $ ivar,n,ifvec,name6)
1142 
1143 c ix is pointer to X, ib is pointer to B
1144 c ixb is pointer to xbar, ibb is pointer to bbar := A*xbar
1145 
1146  call project2_a(x,rvar(ixb,1),rvar(ix,1),rvar(ib,1)
1147  $ ,n,m,mmx,h1,h2,msk,w,ifwt,ifvec,name6)
1148 
1149  ivar(2) = m ! Update number of saved vectors
1150 
1151  tproj = tproj + dnekclock() - etime0
1152 
1153  return
1154  end
1155 c-----------------------------------------------------------------------
1156  subroutine project2_a(x,xbar,xx,bb,n,m,mmx,h1,h2,msk,w,ifwt,ifvec,
1157  $ name6)
1158 
1159  include 'SIZE'
1160 
1161  real x(n),xbar(n),xx(n,1),bb(n,1),h1(n),h2(n),w(n),msk(n)
1162  character*6 name6
1163  logical ifwt,ifvec
1164 
1165  nn = n
1166  if (ifvec) nn=ldim*n
1167 
1168  if (m.gt.0) call add2(x,xbar,n) ! Restore desired solution
1169 
1170  !Uncomment this if using full reorthogonalization
1171 c if (m.eq.mmx) then ! Push old vector off the stack
1172 c do k=2,mmx
1173 c call copy (xx(1,k-1),xx(1,k),nn)
1174 c call copy (bb(1,k-1),bb(1,k),nn)
1175 c enddo
1176 c endif
1177 
1178  m = min(m+1,mmx)
1179  !print *, "m", m
1180  call copy (xx(1,m),x,nn) ! Update (X,B)
1181  call proj_matvec (bb(1,m),xx(1,m),n,h1,h2,msk,name6)
1182  call proj_ortho (xx,bb,n,m,w,ifwt,ifvec,name6) !Update orthogonalization
1183  !Uncomment the if block above if using full reorthogonalization
1184 c call proj_ortho_full (xx,bb,n,m,w,ifwt,ifvec,name6) !Fully reorthogonalize
1185 
1186  return
1187  end
1188 c-----------------------------------------------------------------------
1189  subroutine proj_get_ivar
1190  $ (m,mmx,ixb,ibb,ix,ib,ih1,ih2,ivar,n,ifvec,name6)
1191 
1192  include 'SIZE'
1193  include 'TSTEP'
1194 
1195  logical ifvec
1196  character*6 name6
1197 
1198  integer ivar(10)
1199 
1200  integer icalld
1201  save icalld
1202  data icalld/0/
1203 
1204  m = ivar(2)
1205  mmx = (mxprev-4)/2 ! ivar=0 --> mxprev array
1206  ivar(1) = mmx
1207 
1208  nn = n
1209  if (ifvec) nn = n*ldim ! Number of entries in a vector
1210 
1211 
1212  ih1 = 1
1213  ih2 = ih1 + n
1214  ixb = ih2 + n ! pointer to xbar
1215  ibb = ixb + nn ! " to bbar
1216  ix = ibb + nn ! " to X
1217  ib = ix + nn*mmx ! " to B
1218 
1219  return
1220  end
1221 c-----------------------------------------------------------------------
1222  subroutine laplacep(name,u,mask,mult,ifld,tol,maxi,approx,napprox)
1223 c
1224 c Solve Laplace's equation, with projection onto previous solutions.
1225 c
1226 c Boundary condition strategy:
1227 c
1228 c u = u0 + ub
1229 c
1230 c u0 = 0 on Dirichlet boundaries
1231 c ub = u on Dirichlet boundaries
1232 c
1233 c _
1234 c A ( u0 + ub ) = 0
1235 c
1236 c _ _
1237 c A u0 = - A ub
1238 c
1239 c _ _
1240 c MAM u0 = -M A ub, M is the mask
1241 c
1242 c _
1243 c A u0 = -M A ub , Helmholtz solve with SPD matrix A
1244 c
1245 c u = u0+ub
1246 c
1247  include 'SIZE'
1248  include 'TOTAL'
1249  include 'CTIMER'
1250 c
1251  character*4 name
1252  real u(1),mask(1),mult(1),approx (1)
1253  integer napprox(1)
1254 
1255  parameter(lt=lx1*ly1*lz1*lelt)
1256  common /scrvh/ h1(lt),h2(lt)
1257  common /scruz/ r(lt),ub(lt)
1258 
1259  logical ifstdh
1260  character*4 cname
1261  character*6 name6
1262 
1263  logical ifwt,ifvec
1264 
1265  call chcopy(cname,name,4)
1266  call capit (cname,4)
1267 
1268  call blank (name6,6)
1269  call chcopy(name6,name,4)
1270  ifwt = .true.
1271  ifvec = .false.
1272  isd = 1
1273  imsh = 1
1274  nel = nelfld(ifld)
1275 
1276  n = lx1*ly1*lz1*nel
1277 
1278  call copy (ub,u,n) ! ub = u on boundary
1279  call dsavg(ub) ! Make certain ub is in H1
1280  call rone (h1,n)
1281  call rzero(h2,n)
1282  ! _
1283  call axhelm (r,ub,h1,h2,1,1) ! r = A*ub
1284 
1285  do i=1,n ! _
1286  r(i)=-r(i)*mask(i) ! r = -M*A*ub
1287  enddo
1288 
1289  call dssum (r,lx1,ly1,lz1) ! dssum rhs
1290 
1291  call project1
1292  $ (r,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
1293 
1294  if (nel.eq.nelv) then
1295  call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,binvm1)
1296  else
1297  call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,bintm1)
1298  endif
1299 
1300  call project2
1301  $ (u,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
1302 
1303  call add2(u,ub,n)
1304 
1305  return
1306  end
1307 c-----------------------------------------------------------------------
1308  subroutine givens_rotation(a, b, c, s, r)
1309 
1310  real a, b, c, s, r, h, d
1311 
1312  if(b.ne.0.0) then
1313  h = hypot(a,b) !Can use library version or below implementation
1314  d = 1.0/h
1315  c = abs(a)*d
1316  s = sign(d,a)*b
1317  r = sign(1.0,a)*h
1318  else
1319  c = 1.0
1320  s = 0.0
1321  r = a
1322  endif
1323 
1324  return
1325  end
1326 c-----------------------------------------------------------------------
1327  ! Compilers with Fortran 2008 support should have a
1328  ! library implementation of this (gfortran and ifort do).
1329 
1330  real function hypot(a, b) !Does not handle a = b = 0 case
1331 
1332  real a, b, t, x, c, d, ix
1333 
1334  c = abs(a)
1335  d = abs(b)
1336  x = max(c,d)
1337  ix = 1.0/x
1338  t = min(c,d)
1339  t = ix*t
1340 
1341  hypot = x*sqrt(1.0+t*t)
1342 
1343  return
1344  end
1345 c-----------------------------------------------------------------------
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
chktcg1
subroutine chktcg1(tol, res, h1, h2, mask, mult, imesh, isd)
Definition: hmholtz.f:528
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
cggo
subroutine cggo(x, f, h1, h2, mask, mult, imsh, tin, maxit, isd, binv, name)
Definition: hmholtz.f:612
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
cggo_dg
subroutine cggo_dg(x, f, h1, h2, binv, mask, name, tin, maxit)
Definition: hmholtz.f:1494
dsavg
subroutine dsavg(u)
Definition: ic.f:1878
capit
subroutine capit(lettrs, n)
Definition: ic.f:1648
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
glsc2
function glsc2(x, y, n)
Definition: math.f:794
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
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
copy
subroutine copy(a, b, n)
Definition: math.f:260
blank
subroutine blank(A, N)
Definition: math.f:19
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
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
rone
subroutine rone(a, n)
Definition: math.f:230
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
econj
subroutine econj(kprev, h1, h2, h2inv, ierr)
Definition: navier4.f:201
proj_get_ivar
subroutine proj_get_ivar(m, mmx, ixb, ibb, ix, ib, ih1, ih2, ivar, n, ifvec, name6)
Definition: navier4.f:1191
project1
subroutine project1(b, n, rvar, ivar, h1, h2, msk, w, ifwt, ifvec, name6)
Definition: navier4.f:637
givens_rotation
subroutine givens_rotation(a, b, c, s, r)
Definition: navier4.f:1309
project1_a
subroutine project1_a(xbar, bbar, b, xx, bb, n, m, w, ifwt, ifvec)
Definition: navier4.f:736
hsolve
subroutine hsolve(name, u, r, h1, h2, vmk, vml, imsh, tol, maxit, isd, approx, napprox, bi)
Definition: navier4.f:564
laplacep
subroutine laplacep(name, u, mask, mult, ifld, tol, maxi, approx, napprox)
Definition: navier4.f:1223
updtset
subroutine updtset(p, h1, h2, h2inv, IERR)
Definition: navier4.f:148
chkptol
subroutine chkptol
Definition: navier4.f:269
hypot
real function hypot(a, b)
Definition: navier4.f:1331
proj_ortho
subroutine proj_ortho(xx, bb, n, m, w, ifwt, ifvec, name6)
Definition: navier4.f:853
vlsc3
function vlsc3(X, Y, B, N)
Definition: navier4.f:324
proj_ortho_full_cgs2
subroutine proj_ortho_full_cgs2(xx, bb, n, m, w, ifwt, ifvec, name6)
Definition: navier4.f:1061
hconj
subroutine hconj(approx, k, h1, h2, vml, vmk, ws, name4, ierr)
Definition: navier4.f:434
iproj_chk
function iproj_chk(h1old, h2old, h1, h2, n)
Definition: navier4.f:792
proj_ortho_full_mgs
subroutine proj_ortho_full_mgs(xx, bb, n, m, w, ifwt, ifvec, name6)
Definition: navier4.f:976
project2_a
subroutine project2_a(x, xbar, xx, bb, n, m, mmx, h1, h2, msk, w, ifwt, ifvec, name6)
Definition: navier4.f:1158
setrhs
subroutine setrhs(p, h1, h2, h2inv)
Definition: navier4.f:10
proj_ortho_full
subroutine proj_ortho_full(xx, bb, n, m, w, ifwt, ifvec, name6)
Definition: navier4.f:960
hmhzpf
subroutine hmhzpf(name, u, r, h1, h2, mask, mult, imesh, tli, maxit, isd, bi)
Definition: navier4.f:514
project2
subroutine project2(x, n, rvar, ivar, h1, h2, msk, w, ifwt, ifvec, name6)
Definition: navier4.f:1130
gensoln
subroutine gensoln(p, h1, h2, h2inv)
Definition: navier4.f:111
updrhse
subroutine updrhse(p, h1, h2, h2inv, ierr)
Definition: navier4.f:357
proj_matvec
subroutine proj_matvec(b, x, n, h1, h2, msk, name6)
Definition: navier4.f:828
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341