KTH framework for Nek5000 toolboxes; testing version  0.0.1
subs1.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,binv,
3  $ vol,tin,maxit,matmod)
4 
5 C Conjugate gradient iteration for solution of coupled
6 C Helmholtz equations
7 
8  include 'SIZE'
9  include 'TOTAL'
10  include 'DOMAIN'
11  include 'FDMH1'
12 
13  common /screv/ dpc(lx1*ly1*lz1*lelt)
14  $ , p1(lx1*ly1*lz1*lelt)
15  common /scrch/ p2(lx1*ly1*lz1*lelt)
16  $ , p3(lx1*ly1*lz1*lelt)
17  common /scrsl/ qq1(lx1*ly1*lz1*lelt)
18  $ , qq2(lx1*ly1*lz1*lelt)
19  $ , qq3(lx1*ly1*lz1*lelt)
20  common /scrmg/ pp1(lx1*ly1*lz1*lelt)
21  $ , pp2(lx1*ly1*lz1*lelt)
22  $ , pp3(lx1*ly1*lz1*lelt)
23  $ , wa(lx1*ly1*lz1*lelt)
24  real ap1(1),ap2(1),ap3(1)
25  equivalence(ap1,pp1),(ap2,pp2),(ap3,pp3)
26 
27  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
28  common /cprint/ ifprint
29  logical ifdfrm, iffast, ifh2, ifsolv, ifprint
30 
31  real u1(1),u2(1),u3(1),
32  $ r1(1),r2(1),r3(1),h1(1),h2(1),rmult(1),binv(1)
33 
34 
35  logical iffdm,ifcrsl
36 
37  iffdm = .true.
38  iffdm = .false.
39 c ifcrsl = .true.
40  ifcrsl = .false.
41 
42  nel = nelfld(ifield)
43  nxyz = lx1*ly1*lz1
44  n = nxyz*nel
45 
46  if (istep.le.1.and.iffdm) call set_fdm_prec_h1a
47 
48  tol = tin
49 
50 c overrule input tolerance
51  if (restol(ifield).ne.0) tol=restol(ifield)
52 
53  if (ifcrsl) call set_up_h1_crs_strs(h1,h2,ifield,matmod)
54 
55 c if (nio.eq.0.and.istep.eq.1) write(6,6) ifield,tol,tin
56 c 6 format(i3,1p2e12.4,' ifield, tol, tol_in')
57 
58  if ( .not.ifsolv ) then ! Set logical flags
59  call setfast (h1,h2,imesh)
60  ifsolv = .true.
61  endif
62 
63  call opdot (wa,r1,r2,r3,r1,r2,r3,n)
64  rbnorm = glsc3(wa,binv,rmult,n)
65  rbnorm = sqrt( rbnorm / vol )
66  if (rbnorm .lt. tol**2) then
67  iter = 0
68  r0 = rbnorm
69 c if ( .not.ifprint ) goto 9999
70  if (matmod.ge.0.and.nio.eq.0) write (6,3000)
71  $ istep,iter,rbnorm,r0,tol
72  if (matmod.lt.0.and.nio.eq.0) write (6,3010)
73  $ istep,iter,rbnorm,r0,tol
74  goto 9999
75  endif
76 
77 C Evaluate diagional pre-conidtioner for fluid solve
78  call setprec (dpc,h1,h2,imesh,1)
79  call setprec (wa ,h1,h2,imesh,2)
80  call add2 (dpc,wa,n)
81  if (ldim.eq.3) then
82  call setprec (wa,h1,h2,imesh,3)
83  call add2 (dpc,wa,n)
84  endif
85 c call rone (dpc,n)
86 c call copy (dpc,binv,n)
87 
88  if (iffdm) then
89  call set_fdm_prec_h1b(dpc,h1,h2,nel)
90  call fdm_h1a (pp1,r1,dpc,nel,ktype(1,1,1),wa)
91  call fdm_h1a (pp2,r2,dpc,nel,ktype(1,1,2),wa)
92  call fdm_h1a (pp3,r3,dpc,nel,ktype(1,1,3),wa)
93  call rmask (pp1,pp2,pp3,nel)
94  call opdssum (pp1,pp2,pp3)
95  else
96  call col3 (pp1,dpc,r1,n)
97  call col3 (pp2,dpc,r2,n)
98  if (if3d) call col3 (pp3,dpc,r3,n)
99  endif
100  if (ifcrsl) then
101  call crs_strs(p1,p2,p3,r1,r2,r3)
102  call rmask (p1,p2,p3,nel)
103  else
104  call opzero(p1,p2,p3)
105  endif
106  call opadd2 (p1,p2,p3,pp1,pp2,pp3)
107  rpp1 = op_glsc2_wt(p1,p2,p3,r1,r2,r3,rmult)
108 
109  maxit=200
110  do 1000 iter=1,maxit
111  call axhmsf (ap1,ap2,ap3,p1,p2,p3,h1,h2,matmod)
112  call rmask (ap1,ap2,ap3,nel)
113  call opdssum (ap1,ap2,ap3)
114  pap = op_glsc2_wt(p1,p2,p3,ap1,ap2,ap3,rmult)
115  alpha = rpp1 / pap
116 
117  call opadds (u1,u2,u3,p1 ,p2 ,p3 , alpha,n,2)
118  call opadds (r1,r2,r3,ap1,ap2,ap3,-alpha,n,2)
119 
120  call opdot (wa,r1,r2,r3,r1,r2,r3,n)
121  rbnorm = glsc3(wa,binv,rmult,n)
122  rbnorm = sqrt(rbnorm/vol)
123 
124  if (iter.eq.1) r0 = rbnorm
125 
126  if (rbnorm.lt.tol) then
127  ifin = iter
128  if (nio.eq.0) then
129  if (matmod.ge.0) write(6,3000) istep,ifin,rbnorm,r0,tol
130  if (matmod.lt.0) write(6,3010) istep,ifin,rbnorm,r0,tol
131  endif
132  goto 9999
133  endif
134 
135  if (iffdm) then
136  call fdm_h1a (pp1,r1,dpc,nel,ktype(1,1,1),wa)
137  call fdm_h1a (pp2,r2,dpc,nel,ktype(1,1,2),wa)
138  call fdm_h1a (pp3,r3,dpc,nel,ktype(1,1,3),wa)
139  call rmask (pp1,pp2,pp3,nel)
140  call opdssum (pp1,pp2,pp3)
141  else
142  call col3 (pp1,dpc,r1,n)
143  call col3 (pp2,dpc,r2,n)
144  if (if3d) call col3 (pp3,dpc,r3,n)
145  endif
146 
147  if (ifcrsl) then
148  call crs_strs(qq1,qq2,qq3,r1,r2,r3)
149  call rmask (qq1,qq2,qq3,nel)
150  call opadd2 (pp1,pp2,pp3,qq1,qq2,qq3)
151  endif
152 
153  call opdot (wa,r1,r2,r3,pp1,pp2,pp3,n)
154 
155  rpp2 = rpp1
156  rpp1 = glsc2(wa,rmult,n)
157  beta = rpp1/rpp2
158  call opadds (p1,p2,p3,pp1,pp2,pp3,beta,n,1)
159 
160  1000 continue
161  if (matmod.ge.0.and.nio.eq.0) write (6,3001)
162  $ istep,iter,rbnorm,r0,tol
163  if (matmod.lt.0.and.nio.eq.0) write (6,3011)
164  $ istep,iter,rbnorm,r0,tol
165 
166  9999 continue
167  ifsolv = .false.
168 
169 
170  3000 format(i11,' Helmh3 fluid ',i6,1p3e13.4)
171  3010 format(i11,' Helmh3 mesh ',i6,1p3e13.4)
172  3001 format(i11,' Helmh3 fluid unconverged! ',i6,1p3e13.4)
173  3011 format(i11,' Helmh3 mesh unconverged! ',i6,1p3e13.4)
174 
175  return
176  end
177 c-----------------------------------------------------------------------
178  subroutine setdt
179 c
180 c Set the new time step. All cases covered.
181 c
182  include 'SIZE'
183  include 'SOLN'
184  include 'MVGEOM'
185  include 'INPUT'
186  include 'TSTEP'
187  include 'PARALLEL'
188 
189  common /scruz/ cx(lx1*ly1*lz1*lelt)
190  $ , cy(lx1,ly1,lz1,lelt)
191  $ , cz(lx1,ly1,lz1,lelt)
192 
193  common /cprint/ ifprint
194  logical ifprint
195  common /udxmax/ umax
196  REAL DTOLD
197  SAVE dtold
198  DATA dtold /0.0/
199  REAL DTOpf
200  SAVE dtopf
201  DATA dtopf /0.0/
202  logical iffxdt
203  save iffxdt
204  data iffxdt /.false./
205 C
206 
207  if (param(12).lt.0.or.iffxdt) then
208  iffxdt = .true.
209  param(12) = abs(param(12))
210  dt = param(12)
211  dtopf = dt
212  if (ifmvbd) then
213  call opsub3 (cx,cy,cz,vx,vy,vz,wx,wy,wz)
214  call compute_cfl(umax,cx,cy,cz,1.0)
215  else
216  call compute_cfl(umax,vx,vy,vz,1.0)
217  endif
218 
219  goto 200
220  else IF (param(84).NE.0.0) THEN
221  if (dtold.eq.0.0) then
222  dt =param(84)
223  dtold=param(84)
224  dtopf=param(84)
225  return
226  else
227  dtold=dt
228  dtopf=dt
229  dt=dtopf*param(85)
230  dt=min(dt,param(12))
231  endif
232  endif
233 
234 C
235 C Find DT=DTCFL based on CFL-condition (if applicable)
236 C
237  CALL setdtc
238  dtcfl = dt
239 C
240 C Find DTFS based on surface tension (if applicable)
241 C
242  CALL setdtfs (dtfs)
243 C
244 C Select appropriate DT
245 C
246  IF ((dt.EQ.0.).AND.(dtfs.GT.0.)) THEN
247  dt = dtfs
248  ELSEIF ((dt.GT.0.).AND.(dtfs.GT.0.)) THEN
249  dt = min(dt,dtfs)
250  ELSEIF ((dt.EQ.0.).AND.(dtfs.EQ.0.)) THEN
251  dt = 0.
252  IF (ifflow.AND.nid.EQ.0.AND.ifprint) THEN
253  WRITE (6,*) 'WARNING: CFL-condition & surface tension'
254  WRITE (6,*) ' are not applicable'
255  endif
256  ELSEIF ((dt.GT.0.).AND.(dtfs.EQ.0.)) THEN
257  dt = dt
258  ELSE
259  dt = 0.
260  IF (nio.EQ.0) WRITE (6,*) 'WARNING: DT<0 or DTFS<0'
261  IF (nio.EQ.0) WRITE (6,*) ' Reset DT '
262  endif
263 C
264 C Check DT against user-specified input, DTINIT=PARAM(12).
265 C
266  IF ((dt.GT.0.).AND.(dtinit.GT.0.)) THEN
267  dt = min(dt,dtinit)
268  ELSEIF ((dt.EQ.0.).AND.(dtinit.GT.0.)) THEN
269  dt = dtinit
270  ELSEIF ((dt.GT.0.).AND.(dtinit.EQ.0.)) THEN
271  dt = dt
272  ELSEIF (.not.iffxdt) THEN
273  dt = 0.001
274  IF(nio.EQ.0)WRITE (6,*) 'WARNING: Set DT=0.001 (arbitrarily)'
275  endif
276 C
277 C Check if final time (user specified) has been reached.
278 C
279  200 IF (time+dt .GE. fintim .AND. fintim.NE.0.0) THEN
280 C Last step
281  lastep = 1
282  dt = fintim-time
283  IF (nio.EQ.0) WRITE (6,*) 'Final time step = ',dt
284  endif
285 C
286  courno = dt*umax
287  IF (nio.EQ.0.AND.ifprint.AND.dt.NE.dtold)
288  $ WRITE (6,100) dt,dtcfl,dtfs,dtinit
289  100 FORMAT(5x,'DT/DTCFL/DTFS/DTINIT',4e12.3)
290 C
291 C Put limits on how much DT can change.
292 C
293  IF (dtold.NE.0.0 .AND. lastep.NE.1) THEN
294  dtmin=0.8*dtold
295  dtmax=1.2*dtold
296  dt = min(dtmax,dt)
297  dt = max(dtmin,dt)
298  endif
299  dtold=dt
300 
301 C IF (PARAM(84).NE.0.0) THEN
302 C dt=dtopf*param(85)
303 C dt=min(dt,param(12))
304 C endif
305 
306  if (iffxdt) dt=dtopf
307  courno = dt*umax
308 
309 ! synchronize time step for multiple sessions
310  if (ifneknek) dt = glmin_ms(dt,1)
311 c
312  if (iffxdt.and.abs(courno).gt.10.*abs(ctarg)) then
313  if (nid.eq.0) write(6,*) 'CFL, Ctarg!',courno,ctarg
314  call emerxit
315  endif
316 
317 
318  return
319  end
320 C
321 C--------------------------------------------------------
322 C
323  subroutine cvgnlps (ifconv)
324 C----------------------------------------------------------------------
325 C
326 C Check convergence for non-linear passisve scalar solver.
327 C Relevant for solving heat transport problems with radiation b.c.
328 C
329 C----------------------------------------------------------------------
330  include 'SIZE'
331  include 'INPUT'
332  include 'TSTEP'
333  LOGICAL IFCONV
334 C
335  IF (ifnonl(ifield)) THEN
336  ifconv = .false.
337  ELSE
338  ifconv = .true.
339  return
340  endif
341 C
342  tnorm1 = tnrmh1(ifield-1)
343  CALL unorm
344  tnorm2 = tnrmh1(ifield-1)
345  eps = abs((tnorm2-tnorm1)/tnorm2)
346  IF (eps .LT. tolnl) ifconv = .true.
347 C
348  return
349  end
350 C
351  subroutine unorm
352 C---------------------------------------------------------------------
353 C
354 C Norm calculation.
355 C
356 C---------------------------------------------------------------------
357  include 'SIZE'
358  include 'SOLN'
359  include 'TSTEP'
360 C
361  IF (ifield.EQ.1) THEN
362 C
363 C Compute norms of the velocity.
364 C Compute time mean (L2) of the inverse of the time step.
365 C Compute L2 in time, H1 in space of the velocity.
366 C
367  CALL normvc (vnrmh1,vnrmsm,vnrml2,vnrml8,vx,vy,vz)
368  IF (istep.EQ.0) return
369  IF (istep.EQ.1) THEN
370  dtinvm = 1./dt
371  vmean = vnrml8
372  ELSE
373  tden = time
374  if (time.le.0) tden = abs(time)+1.e-9
375  arg = ((time-dt)*dtinvm**2+1./dt)/tden
376  if (arg.gt.0) dtinvm = sqrt(arg)
377  arg = ((time-dt)*vmean**2+dt*vnrmh1**2)/tden
378  if (arg.gt.0) vmean = sqrt(arg)
379  endif
380  ELSE
381 C
382 C Compute norms of a passive scalar
383 C
384  CALL normsc (tnrmh1(ifield-1),tnrmsm(ifield-1),
385  $ tnrml2(ifield-1),tnrml8(ifield-1),
386  $ t(1,1,1,1,ifield-1),imesh)
387  tmean(ifield-1) = 0.
388  endif
389 C
390  return
391  end
392 C
393  subroutine chktmg (tol,res,w1,w2,mult,mask,imesh)
394 C-------------------------------------------------------------------
395 C
396 C Check that the tolerances are not too small for the MG-solver.
397 C Important when calling the MG-solver (Gauss-Lobatto mesh).
398 C Note: direct stiffness summation
399 C
400 C-------------------------------------------------------------------
401  include 'SIZE'
402  include 'INPUT'
403  include 'MASS'
404  include 'EIGEN'
405 C
406  REAL RES (LX1,LY1,LZ1,1)
407  REAL W1 (LX1,LY1,LZ1,1)
408  REAL W2 (LX1,LY1,LZ1,1)
409  REAL MULT (LX1,LY1,LZ1,1)
410  REAL MASK (LX1,LY1,LZ1,1)
411 C
412 C Single or double precision???
413 C
414  delta = 1.e-9
415  x = 1.+delta
416  y = 1.
417  diff = abs(x-y)
418  IF (diff.EQ.0.) eps = 1.e-6*eigga/eigaa
419  IF (diff.GT.0.) eps = 1.e-13*eigga/eigaa
420 C
421  IF (imesh.EQ.1) nl = nelv
422  IF (imesh.EQ.2) nl = nelt
423  ntot1 = lx1*ly1*lz1*nl
424  CALL copy (w1,res,ntot1)
425 C
426  CALL dssum (w1,lx1,ly1,lz1)
427 C
428  IF (imesh.EQ.1) THEN
429  CALL col3 (w2,binvm1,w1,ntot1)
430  rinit = sqrt(glsc3(w2,w1,mult,ntot1)/volvm1)
431  ELSE
432  CALL col3 (w2,bintm1,w1,ntot1)
433  rinit = sqrt(glsc3(w2,w1,mult,ntot1)/voltm1)
434  endif
435  rmin = eps*rinit
436  IF (tol.LT.rmin) THEN
437  tolold=tol
438  tol = rmin
439  IF (nio.EQ.0)
440  $ WRITE (6,*) 'New MG-tolerance (RINIT*epsm*cond) = ',tol,tolold
441  endif
442 C
443  CALL rone (w1,ntot1)
444  bcneu1 = glsc3(w1,mask,mult,ntot1)
445  bcneu2 = glsc3(w1,w1 ,mult,ntot1)
446  bctest = abs(bcneu1-bcneu2)
447  IF (bctest .LT. .1) THEN
448  otr = glsum(res,ntot1)
449  tolmin = abs(otr)*eigga/eigaa
450  IF (tol .LT. tolmin) THEN
451  tolold = tol
452  tol = tolmin
453  IF (nio.EQ.0)
454  $ WRITE (6,*) 'New MG-tolerance (OTR) = ',tol,tolold
455  endif
456  endif
457 C
458  return
459  end
460 C
461 C
462  subroutine setdtc
463 C--------------------------------------------------------------
464 C
465 C Compute new timestep based on CFL-condition
466 C
467 C--------------------------------------------------------------
468  include 'SIZE'
469  include 'GEOM'
470  include 'MVGEOM'
471  include 'MASS'
472  include 'INPUT'
473  include 'SOLN'
474  include 'TSTEP'
475 C
476  common /ctmp1/ u(lx1,ly1,lz1,lelv)
477  $ , v(lx1,ly1,lz1,lelv)
478  $ , w(lx1,ly1,lz1,lelv)
479  common /ctmp0/ x(lx1,ly1,lz1,lelv)
480  $ , r(lx1,ly1,lz1,lelv)
481  common /udxmax/ umax
482 
483  common /scruz/ cx(lx1*ly1*lz1*lelv)
484  $ , cy(lx1,ly1,lz1,lelv)
485  $ , cz(lx1,ly1,lz1,lelv)
486 C
487 C
488  REAL VCOUR
489  SAVE vcour
490 C
491  INTEGER IFIRST
492  SAVE ifirst
493  DATA ifirst/0/
494 C
495 C
496 C Steady state => all done
497 C
498 C
499  IF (.NOT.iftran) THEN
500  ifirst=1
501  lastep=1
502  return
503  endif
504 C
505  irst = param(46)
506  if (irst.gt.0) ifirst=1
507 C
508 C First time around
509 C
510  IF (ifirst.EQ.0) THEN
511  dt = dtinit
512  IF (ifflow) THEN
513  ifield = 1
514  CALL bcdirvc (vx,vy,vz,v1mask,v2mask,v3mask)
515  endif
516  endif
517  ifirst=ifirst+1
518 C
519  dtold = dt
520 C
521 C Convection ?
522 C
523 C Don't enforce Courant condition if there is no convection.
524 C
525 C
526  iconv=0
527  IF (ifflow .AND. ifnav) iconv=1
528  IF (ifwcno) iconv=1
529  IF (ifheat) THEN
530  DO 10 ipscal=0,npscal
531  IF (ifadvc(ipscal+2)) iconv=1
532  10 CONTINUE
533  endif
534 
535  IF (iconv.EQ.0) THEN
536  dt=0.
537  return
538  endif
539 C
540 C
541 C Find Courant and Umax
542 C
543 C
544  ntot = lx1*ly1*lz1*nelv
545  ntotl = lx1*ly1*lz1*lelv
546  ntotd = ntotl*ldim
547  cold = courno
548  cmax = 1.2*ctarg
549  cmin = 0.8*ctarg
550 
551  if (ifmvbd) then
552  call opsub3 (cx,cy,cz,vx,vy,vz,wx,wy,wz)
553  call compute_cfl(umax,cx,cy,cz,1.0)
554  else
555  call compute_cfl(umax,vx,vy,vz,1.0)
556  endif
557 
558 c if (nio.eq.0) write(6,1) istep,time,umax,cmax
559 c 1 format(i9,1p3e12.4,' cmax')
560 
561 C Zero DT
562 
563  IF (dt .EQ. 0.0) THEN
564 
565  IF (umax .NE. 0.0) THEN
566  dt = ctarg/umax
567  vcour = umax
568  ELSEIF (ifflow) THEN
569 C
570 C We'll use the body force to predict max velocity
571 C
572  CALL setprop
573  ifield = 1
574 C
575  CALL makeuf
576  CALL opdssum (bfx,bfy,bfz)
577  CALL opcolv (bfx,bfy,bfz,binvm1)
578  fmax=0.0
579  CALL rzero (u,ntotd)
580  DO 600 i=1,ntot
581  u(i,1,1,1) = abs(bfx(i,1,1,1))
582  v(i,1,1,1) = abs(bfy(i,1,1,1))
583  w(i,1,1,1) = abs(bfz(i,1,1,1))
584  600 CONTINUE
585  fmax = glmax(u,ntotd)
586  density = avtran(1)
587  amax = fmax/density
588  dxchar = sqrt( (xm1(1,1,1,1)-xm1(2,1,1,1))**2 +
589  $ (ym1(1,1,1,1)-ym1(2,1,1,1))**2 +
590  $ (zm1(1,1,1,1)-zm1(2,1,1,1))**2 )
591  dxchar = glmin(dxchar,1)
592  IF (amax.NE.0.) THEN
593  dt = sqrt(ctarg*dxchar/amax)
594  ELSE
595  IF (nid.EQ.0)
596  $ WRITE (6,*) 'CFL: Zero velocity and body force'
597  dt = 0.0
598  return
599  endif
600  ELSEIF (ifwcno) THEN
601  IF (nid.EQ.0)
602  $ WRITE (6,*) ' Stefan problem with no fluid flow'
603  dt = 0.0
604  return
605  endif
606 C
607  ELSEIF ((dt.GT.0.0).AND.(umax.NE.0.0)) THEN
608 C
609 C
610 C Nonzero DT & nonzero velocity
611 C
612 C
613  courno = dt*umax
614  vold = vcour
615  vcour = umax
616  IF (ifirst.EQ.1) THEN
617  cold = courno
618  vold = vcour
619  endif
620  cpred = 2.*courno-cold
621 C
622 C Change DT if it is too big or if it is too small
623 C
624 c if (nid.eq.0)
625 c $write(6,917) dt,umax,vold,vcour,cpred,cmax,courno,cmin
626 c 917 format(' dt',4f9.5,4f10.6)
627  IF(courno.GT.cmax .OR. cpred.GT.cmax .OR. courno.LT.cmin) THEN
628 C
629  a=(vcour-vold)/dt
630  b=vcour
631 C -C IS Target Courant number
632  c=-ctarg
633  discr=b**2-4*a*c
634  dtold=dt
635  IF(discr.LE.0.0)THEN
636 c if (nio.eq.0)
637 c $ PRINT*,'Problem calculating new DT Discriminant=',discr
638  dt=dt*(ctarg/courno)
639 C IF(DT.GT.DTOLD) DT=DTOLD
640  ELSE IF(abs((vcour-vold)/vcour).LT.0.001)THEN
641 C Easy: same v as before (LINEARIZED)
642  dt=dt*(ctarg/courno)
643 c if (nid.eq.0)
644 c $write(6,918) dt,dthi,dtlow,discr,a,b,c
645 c 918 format(' d2',4f9.5,4f10.6)
646  ELSE
647  dtlow=(-b+sqrt(discr) )/(2.0*a)
648  dthi =(-b-sqrt(discr) )/(2.0*a)
649  IF(dthi .GT. 0.0 .AND. dtlow .GT. 0.0)THEN
650  dt = min(dthi,dtlow)
651 c if (nid.eq.0)
652 c $write(6,919) dt,dthi,dtlow,discr,a,b,c
653 c 919 format(' d3',4f9.5,4f10.6)
654  ELSE IF(dthi .LE. 0.0 .AND. dtlow .LE. 0.0)THEN
655 c PRINT*,'DTLOW,DTHI',DTLOW,DTHI
656 c PRINT*,'WARNING: Abnormal DT from CFL-condition'
657 c PRINT*,' Keep going'
658  dt=dt*(ctarg/courno)
659  ELSE
660 C Normal case; 1 positive root, one negative root
661  dt = max(dthi,dtlow)
662 c if (nid.eq.0)
663 c $write(6,929) dt,dthi,dtlow,discr,a,b,c
664 c 929 format(' d4',4f9.5,4f10.6)
665  endif
666  endif
667 C We'll increase gradually-- make it the geometric mean between
668 c if (nid.eq.0)
669 c $write(6,939) dt,dtold
670 c 939 format(' d5',4f9.5,4f10.6)
671  IF (dtold/dt .LT. 0.2) dt = dtold*5
672  endif
673 C
674  endif
675 C
676  return
677  end
678 C
679  subroutine setdtfs (dtfs)
680 C
681  include 'SIZE'
682  include 'INPUT'
683  include 'SOLN'
684  include 'GEOM'
685  include 'TSTEP'
686  common /ctmp0/ stc(lx1,ly1,lz1),sigst(lx1,ly1),dtst(lx1,ly1)
687  character cb*3,cb2*2
688 C
689 C Applicable?
690 C
691  IF (.NOT.ifsurt) THEN
692  dtfs = 0.0
693  return
694  endif
695 C
696  nface = 2*ldim
697  nxz1 = lx1*lz1
698  ntot1 = lx1*ly1*lz1*nelv
699  dtfs = 1.1e+10
700 C
701 C Kludge : for debugging purpose Fudge Factor comes in
702 C as PARAM(45)
703 C
704  factor = param(45)
705  IF (factor .LE. 0.0) factor=1.0
706  factor = factor / sqrt( pi**3 )
707 C
708  IF (istep.EQ.1) CALL setprop
709 C
710  ifield = 1
711  rhomin = glmin( vtrans(1,1,1,1,ifield),ntot1 )
712 C
713  DO 100 iel=1,nelv
714  DO 100 ifc=1,nface
715  cb =cbc(ifc,iel,ifield)
716  cb2=cbc(ifc,iel,ifield)
717  IF (cb2.NE.'MS' .AND. cb2.NE.'ms') GOTO 100
718  IF (cb2(1:1).EQ.'M') THEN
719  bc4 = bc(4,ifc,iel,ifield)
720  CALL cfill (sigst,bc4,nxz1)
721  ELSE
722  CALL faceis (cb,stc,iel,ifc,lx1,ly1,lz1)
723  CALL facexs (sigst,stc,ifc,0)
724  endif
725  sigmax = vlmax(sigst,nxz1)
726  IF (sigmax.LE.0.0) GOTO 100
727  rhosig = sqrt( rhomin/sigmax )
728  IF (ldim.EQ.2) THEN
729  CALL cdxmin2 (dtst,rhosig,iel,ifc,ifaxis)
730  ELSE
731  CALL cdxmin3 (dtst,rhosig,iel,ifc)
732  endif
733  dtmin = vlmin( dtst,nxz1 )
734  dtfs = min( dtfs,dtmin )
735  100 CONTINUE
736  dtfs = glmin( dtfs,1 )
737 C
738  IF (dtfs .GT. 1.e+10) dtfs = 0.0
739  dtfs = dtfs * factor
740 C
741  IF (dtfs.EQ.0.0) THEN
742  IF (istep.EQ.1.AND.nio.EQ.0) THEN
743  WRITE (6,*) ' Warning - zero surface-tension may results in'
744  WRITE (6,*) ' instability of free-surface update'
745  endif
746  endif
747 C
748  return
749  end
750  subroutine cdxmin2 (dtst,rhosig,iel,ifc,ifaxis)
751 C
752  include 'SIZE'
753  include 'GEOM'
754  include 'DXYZ'
755  common /delrst/ drst(lx1),drsti(lx1)
756  common /ctmp0/ xfm1(lx1),yfm1(lx1),t1xf(lx1),t1yf(lx1)
757  dimension dtst(lx1,1)
758  LOGICAL IFAXIS
759 C
760  delta = 1.e-9
761  x = 1.+delta
762  y = 1.
763  diff = abs(x-y)
764  IF (diff.EQ.0.) eps = 1.e-6
765  IF (diff.GT.0.) eps = 1.e-13
766 C
767  CALL facec2 (xfm1,yfm1,xm1(1,1,1,iel),ym1(1,1,1,iel),ifc)
768 C
769  IF (ifc.EQ.1 .OR. ifc.EQ.3) THEN
770  CALL mxm (dxm1,lx1,xfm1,lx1,t1xf,1)
771  CALL mxm (dxm1,lx1,yfm1,lx1,t1yf,1)
772  ELSE
773  IF (ifaxis) CALL setaxdy ( ifrzer(iel) )
774  CALL mxm (dym1,ly1,xfm1,ly1,t1xf,1)
775  CALL mxm (dym1,ly1,yfm1,ly1,t1yf,1)
776  endif
777 C
778  IF (ifaxis) THEN
779  DO 100 ix=1,lx1
780  IF (yfm1(ix) .LT. eps) THEN
781  dtst(ix,1) = 1.e+10
782  ELSE
783  xj = sqrt( t1xf(ix)**2 + t1yf(ix)**2 )*drst(ix)
784  dtst(ix,1) = rhosig * sqrt( xj**3 ) * yfm1(ix)
785  endif
786  100 CONTINUE
787  ELSE
788  DO 200 ix=1,lx1
789  xj = sqrt( t1xf(ix)**2 + t1yf(ix)**2 )*drst(ix)
790  dtst(ix,1) = rhosig * sqrt( xj**3 )
791  200 CONTINUE
792  endif
793 C
794  return
795  end
796  subroutine cdxmin3 (dtst,rhosig,iel,ifc)
797 C
798  include 'SIZE'
799  include 'GEOM'
800  include 'DXYZ'
801  common /delrst/ drst(lx1),drsti(lx1)
802  common /ctmp0/ xfm1(lx1,ly1),yfm1(lx1,ly1),zfm1(lx1,ly1)
803  common /ctmp1/ drm1(lx1,lx1),drtm1(lx1,ly1)
804  $ , dsm1(lx1,lx1),dstm1(lx1,ly1)
805  common /scrmg/ xrm1(lx1,ly1),yrm1(lx1,ly1),zrm1(lx1,ly1)
806  $ , xsm1(lx1,ly1),ysm1(lx1,ly1),zsm1(lx1,ly1)
807  dimension dtst(lx1,ly1)
808 C
809  call facexv (xfm1,yfm1,zfm1,xm1(1,1,1,iel),ym1(1,1,1,iel),
810  $ zm1(1,1,1,iel),ifc,0)
811  call setdrs (drm1,drtm1,dsm1,dstm1,ifc)
812 C
813  CALL mxm (drm1,lx1, xfm1,lx1,xrm1,ly1)
814  CALL mxm (drm1,lx1, yfm1,lx1,yrm1,ly1)
815  CALL mxm (drm1,lx1, zfm1,lx1,zrm1,ly1)
816  CALL mxm (xfm1,lx1,dstm1,ly1,xsm1,ly1)
817  CALL mxm (yfm1,lx1,dstm1,ly1,ysm1,ly1)
818  CALL mxm (zfm1,lx1,dstm1,ly1,zsm1,ly1)
819 C
820  DO 100 ix=1,lx1
821  DO 100 iy=1,ly1
822  delr = xrm1(ix,iy)**2 + yrm1(ix,iy)**2 + zrm1(ix,iy)**2
823  dels = xsm1(ix,iy)**2 + ysm1(ix,iy)**2 + zsm1(ix,iy)**2
824  delr = sqrt( delr )*drst(ix)
825  dels = sqrt( dels )*drst(iy)
826  xj = min( delr,dels )
827  dtst(ix,iy) = rhosig * sqrt( xj**3 )
828  100 CONTINUE
829 C
830  return
831  end
832 C
833  FUNCTION facdot(A,B,IFACE1)
834 C
835 C Take the dot product of A and B on the surface IFACE1 of element IE.
836 C
837 C IFACE1 is in the preprocessor notation
838 C IFACE is the dssum notation.
839 C 5 Jan 1989 15:12:22 PFF
840 C
841  include 'SIZE'
842  include 'TOPOL'
843  dimension a(lx1,ly1,lz1),b(lx1,ly1)
844 C
845 C Set up counters
846 C
847  CALL dsset(lx1,ly1,lz1)
848  iface = eface1(iface1)
849  js1 = skpdat(1,iface)
850  jf1 = skpdat(2,iface)
851  jskip1 = skpdat(3,iface)
852  js2 = skpdat(4,iface)
853  jf2 = skpdat(5,iface)
854  jskip2 = skpdat(6,iface)
855 C
856  sum=0.0
857  i = 0
858  DO 100 j2=js2,jf2,jskip2
859  DO 100 j1=js1,jf1,jskip1
860  i = i+1
861  sum = sum + a(j1,j2,1)*b(i,1)
862  100 CONTINUE
863 C
864  facdot = sum
865 C
866  return
867  end
868 C
869  subroutine fcaver(xaver,a,iel,iface1)
870 C------------------------------------------------------------------------
871 C
872 C Compute the average of A over the face IFACE1 in element IEL.
873 C
874 C A is a (NX,NY,NZ) data structure
875 C IFACE1 is in the preprocessor notation
876 C IFACE is the dssum notation.
877 C------------------------------------------------------------------------
878  include 'SIZE'
879  include 'GEOM'
880  include 'TOPOL'
881  REAL A(LX1,LY1,LZ1,1)
882 C
883  fcarea = 0.
884  xaver = 0.
885 C
886 C Set up counters
887 C
888  CALL dsset(lx1,ly1,lz1)
889  iface = eface1(iface1)
890  js1 = skpdat(1,iface)
891  jf1 = skpdat(2,iface)
892  jskip1 = skpdat(3,iface)
893  js2 = skpdat(4,iface)
894  jf2 = skpdat(5,iface)
895  jskip2 = skpdat(6,iface)
896 C
897  i = 0
898  DO 100 j2=js2,jf2,jskip2
899  DO 100 j1=js1,jf1,jskip1
900  i = i+1
901  fcarea = fcarea+area(i,1,iface1,iel)
902  xaver = xaver +area(i,1,iface1,iel)*a(j1,j2,1,iel)
903  100 CONTINUE
904 C
905  xaver = xaver/fcarea
906  return
907  end
908  subroutine faccl2(a,b,iface1)
909 C
910 C Collocate B with A on the surface IFACE1 of element IE.
911 C
912 C A is a (NX,NY,NZ) data structure
913 C B is a (NX,NY,IFACE) data structure
914 C IFACE1 is in the preprocessor notation
915 C IFACE is the dssum notation.
916 C 5 Jan 1989 15:12:22 PFF
917 C
918  include 'SIZE'
919  include 'TOPOL'
920  dimension a(lx1,ly1,lz1),b(lx1,ly1)
921 C
922 C Set up counters
923 C
924  CALL dsset(lx1,ly1,lz1)
925  iface = eface1(iface1)
926  js1 = skpdat(1,iface)
927  jf1 = skpdat(2,iface)
928  jskip1 = skpdat(3,iface)
929  js2 = skpdat(4,iface)
930  jf2 = skpdat(5,iface)
931  jskip2 = skpdat(6,iface)
932 C
933  i = 0
934  DO 100 j2=js2,jf2,jskip2
935  DO 100 j1=js1,jf1,jskip1
936  i = i+1
937  a(j1,j2,1) = a(j1,j2,1)*b(i,1)
938  100 CONTINUE
939 C
940  return
941  end
942 C
943  subroutine faccl3(a,b,c,iface1)
944 C
945 C Collocate B with A on the surface IFACE1 of element IE.
946 C
947 C A is a (NX,NY,NZ) data structure
948 C B is a (NX,NY,IFACE) data structure
949 C IFACE1 is in the preprocessor notation
950 C IFACE is the dssum notation.
951 C 5 Jan 1989 15:12:22 PFF
952 C
953  include 'SIZE'
954  include 'TOPOL'
955  dimension a(lx1,ly1,lz1),b(lx1,ly1,lz1),c(lx1,ly1)
956 C
957 C Set up counters
958 C
959  CALL dsset(lx1,ly1,lz1)
960  iface = eface1(iface1)
961  js1 = skpdat(1,iface)
962  jf1 = skpdat(2,iface)
963  jskip1 = skpdat(3,iface)
964  js2 = skpdat(4,iface)
965  jf2 = skpdat(5,iface)
966  jskip2 = skpdat(6,iface)
967 C
968  i = 0
969  DO 100 j2=js2,jf2,jskip2
970  DO 100 j1=js1,jf1,jskip1
971  i = i+1
972  a(j1,j2,1) = b(j1,j2,1)*c(i,1)
973  100 CONTINUE
974 C
975  return
976  end
977  subroutine faddcl3(a,b,c,iface1)
978 C
979 C Collocate B with C and add to A on the surface IFACE1 of element IE.
980 C
981 C A is a (NX,NY,NZ) data structure
982 C B is a (NX,NY,NZ) data structure
983 C C is a (NX,NY,IFACE) data structure
984 C IFACE1 is in the preprocessor notation
985 C IFACE is the dssum notation.
986 C 29 Jan 1990 18:00 PST PFF
987 C
988  include 'SIZE'
989  include 'TOPOL'
990  dimension a(lx1,ly1,lz1),b(lx1,ly1,lz1),c(lx1,ly1)
991 C
992 C Set up counters
993 C
994  CALL dsset(lx1,ly1,lz1)
995  iface = eface1(iface1)
996  js1 = skpdat(1,iface)
997  jf1 = skpdat(2,iface)
998  jskip1 = skpdat(3,iface)
999  js2 = skpdat(4,iface)
1000  jf2 = skpdat(5,iface)
1001  jskip2 = skpdat(6,iface)
1002 C
1003  i = 0
1004  DO 100 j2=js2,jf2,jskip2
1005  DO 100 j1=js1,jf1,jskip1
1006  i = i+1
1007  a(j1,j2,1) = a(j1,j2,1) + b(j1,j2,1)*c(i,1)
1008  100 CONTINUE
1009 C
1010  return
1011  end
1012 c-----------------------------------------------------------------------
1013  subroutine sethlm (h1,h2,intloc)
1014 
1015 c Set the variable property arrays H1 and H2
1016 c in the Helmholtz equation.
1017 c (associated with variable IFIELD)
1018 c INTLOC = integration type
1019 
1020  include 'SIZE'
1021  include 'INPUT'
1022  include 'SOLN'
1023  include 'TSTEP'
1024 
1025  real h1(1),h2(1)
1026 
1027  nel = nelfld(ifield)
1028  ntot1 = lx1*ly1*lz1*nel
1029 
1030  if (iftran) then
1031  dtbd = bd(1)/dt
1032  call copy (h1,vdiff(1,1,1,1,ifield),ntot1)
1033  if (intloc.eq.0) then
1034  call rzero (h2,ntot1)
1035  else
1036  call cmult2 (h2,vtrans(1,1,1,1,ifield),dtbd,ntot1)
1037  endif
1038 
1039 c if (ifield.eq.1 .and. ifanls) then ! this should be replaced
1040 c const = 2. ! with a correct stress
1041 c call cmult (h1,const,ntot1) ! formulation
1042 c endif
1043 
1044  ELSE
1045  CALL copy (h1,vdiff(1,1,1,1,ifield),ntot1)
1046  CALL rzero (h2,ntot1)
1047  endif
1048 
1049  return
1050  end
1051 c-----------------------------------------------------------------------
1052  subroutine nekuvp (iel)
1053 C------------------------------------------------------------------
1054 C
1055 C Generate user-specified material properties
1056 C
1057 C------------------------------------------------------------------
1058  include 'SIZE'
1059  include 'INPUT'
1060  include 'SOLN'
1061  include 'TSTEP'
1062  include 'PARALLEL'
1063  include 'NEKUSE'
1064  ielg = lglel(iel)
1065 c IF (IFSTRS .AND. IFIELD.EQ.1) CALL STNRINV ! don't call! pff, 2007
1066  DO 10 k=1,lz1
1067  DO 10 j=1,ly1
1068  DO 10 i=1,lx1
1069  if (optlevel.le.2) CALL nekasgn (i,j,k,iel)
1070  CALL uservp (i,j,k,ielg)
1071  vdiff(i,j,k,iel,ifield) = udiff
1072  vtrans(i,j,k,iel,ifield) = utrans
1073  10 CONTINUE
1074  return
1075  end
1076 C
1077  subroutine diagnos
1078  return
1079  end
1080 
1081 c-----------------------------------------------------------------------
1082  subroutine setsolv
1083  include 'SIZE'
1084 
1085  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1086  logical ifdfrm, iffast, ifh2, ifsolv
1087 
1088  ifsolv = .false.
1089 
1090  return
1091  end
1092 c-----------------------------------------------------------------------
1093  subroutine hmhzsf (name,u1,u2,u3,r1,r2,r3,h1,h2,
1094  $ rmask1,rmask2,rmask3,rmult,
1095  $ tol,maxit,matmod)
1096 
1097 c Solve coupled Helmholtz equations (stress formulation)
1098 
1099 
1100  include 'SIZE'
1101  include 'INPUT'
1102  include 'MASS'
1103  include 'SOLN' ! For outpost diagnostic call
1104  include 'TSTEP'
1105  include 'ORTHOSTRS'
1106  include 'CTIMER'
1107 
1108  real u1(1),u2(1),u3(1),r1(1),r2(1),r3(1),h1(1),h2(1)
1109  real rmask1(1),rmask2(1),rmask3(1),rmult(1)
1110  character name*4
1111 
1112 #ifdef TIMER
1113  nhmhz = nhmhz + 1
1114  etime1 = dnekclock()
1115 #endif
1116 
1117  nel = nelfld(ifield)
1118  vol = volfld(ifield)
1119  n = lx1*ly1*lz1*nel
1120 
1121  napproxstrs(1) = 0
1122  iproj = 0
1123  if (ifprojfld(ifield)) iproj = param(94)
1124  if (iproj.gt.0.and.istep.ge.iproj) napproxstrs(1)=param(93)
1125  napproxstrs(1)=min(napproxstrs(1),mxprev)
1126 
1127  call rmask (r1,r2,r3,nel)
1128  call opdssum (r1,r2,r3)
1129  call rzero3 (u1,u2,u3,n)
1130 
1131  if (imesh.eq.1) then
1132  call chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,binvm1
1133  $ ,vol,tol,nel)
1134 
1135  call strs_project_a(r1,r2,r3,h1,h2,rmult,ifield,ierr,matmod)
1136 
1137  call cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,binvm1
1138  $ ,vol,tol,maxit,matmod)
1139 
1140  call strs_project_b(u1,u2,u3,h1,h2,rmult,ifield,ierr)
1141 
1142  else
1143 
1144  call chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,bintm1
1145  $ ,vol,tol,nel)
1146  call cggosf (u1,u2,u3,r1,r2,r3,h1,h2,rmult,bintm1
1147  $ ,vol,tol,maxit,matmod)
1148 
1149  endif
1150 
1151 #ifdef TIMER
1152  thmhz=thmhz+(dnekclock()-etime1)
1153 #endif
1154 
1155  return
1156  end
1157 
1158  subroutine chktcgs (r1,r2,r3,rmask1,rmask2,rmask3,rmult,binv,
1159  $ vol,tol,nel)
1160 C-------------------------------------------------------------------
1161 C
1162 C Check that the tolerances are not too small for the CG-solver.
1163 C Important when calling the CG-solver (Gauss-Lobatto mesh) with
1164 C zero Neumann b.c.
1165 C
1166 C-------------------------------------------------------------------
1167  include 'SIZE'
1168  include 'INPUT'
1169  include 'MASS'
1170  include 'EIGEN'
1171  common /cprint/ ifprint
1172  logical ifprint
1173  common /ctmp0/ wa(lx1,ly1,lz1,lelt)
1174 C
1175  dimension r1(lx1,ly1,lz1,1)
1176  $ , r2(lx1,ly1,lz1,1)
1177  $ , r3(lx1,ly1,lz1,1)
1178  $ , rmask1(lx1,ly1,lz1,1)
1179  $ , rmask2(lx1,ly1,lz1,1)
1180  $ , rmask3(lx1,ly1,lz1,1)
1181  $ , rmult(lx1,ly1,lz1,1)
1182  $ , binv(lx1,ly1,lz1,1)
1183 C
1184  ntot1 = lx1*ly1*lz1*nel
1185 C
1186  IF (eigaa .NE. 0.0) THEN
1187  acondno = eigga/eigaa
1188  ELSE
1189  acondno = 10.0
1190  endif
1191 C
1192 C Check Single or double precision
1193 C
1194  delta = 1.0e-9
1195  x = 1.0 + delta
1196  y = 1.0
1197  diff = abs(x - y)
1198  IF (diff .EQ. 0.0) eps = 1.0e-6
1199  IF (diff .GT. 0.0) eps = 1.0e-13
1200 C
1201  CALL opdot (wa,r1,r2,r3,r1,r2,r3,ntot1)
1202  rinit = glsc3(wa,binv,rmult,ntot1)
1203  rinit = sqrt(rinit/vol)
1204  rmin = eps*rinit
1205 C
1206  IF (tol.LT.rmin) THEN
1207  tolold = tol
1208  tol = rmin
1209 c IF (NIO.EQ.0 .AND. IFPRINT)
1210 c $ WRITE(6,*)'New CG1(stress)-tolerance (RINIT*epsm) = ',TOL,TOLOLD
1211  endif
1212 C
1213  IF (ldim.EQ.2) THEN
1214  CALL add3 (wa,rmask1,rmask2,ntot1)
1215  ELSE
1216  CALL add4 (wa,rmask1,rmask2,rmask3,ntot1)
1217  endif
1218  bcneu1 = glsc2(wa,rmult,ntot1)
1219  bcneu2 = (ldim) * glsum(rmult,ntot1)
1220  bctest = abs(bcneu1 - bcneu2)
1221  IF (bctest .LT. 0.1) THEN
1222  IF (ldim.EQ.2) THEN
1223  CALL add3 (wa,r1,r2,ntot1)
1224  ELSE
1225  CALL add4 (wa,r1,r2,r3,ntot1)
1226  endif
1227  otr = glsc2(wa,rmult,ntot1) / ( ldim )
1228  tolmin = abs(otr) * acondno
1229  IF (tol .LT. tolmin) THEN
1230  tolold = tol
1231  tol = tolmin
1232 c IF (NIO.EQ.0)
1233 c $ WRITE (6,*) 'New CG1(stress)-tolerance (OTR) = ',TOL,TOLOLD
1234  endif
1235  endif
1236 C
1237  return
1238  end
1239 c-----------------------------------------------------------------------
1240  subroutine axhmsf (au1,au2,au3,u1,u2,u3,h1,h2,matmod)
1241 
1242 C Compute the coupled Helmholtz matrix-vector products
1243 
1244 C Fluid (MATMOD .GE. 0) : Hij Uj = Aij*Uj + H2*B*Ui
1245 C Solid (MATMOD .LT. 0) : Hij Uj = Kij*Uj
1246 C
1247 C-----------------------------------------------------------------------
1248  include 'SIZE'
1249  include 'INPUT'
1250  include 'GEOM'
1251  include 'MASS'
1252  include 'TSTEP'
1253  include 'CTIMER'
1254 
1255  common /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
1256  logical ifdfrm, iffast, ifh2, ifsolv
1257 C
1258  dimension au1(lx1,ly1,lz1,1)
1259  $ , au2(lx1,ly1,lz1,1)
1260  $ , au3(lx1,ly1,lz1,1)
1261  $ , u1(lx1*ly1*lz1*1)
1262  $ , u2(lx1*ly1*lz1*1)
1263  $ , u3(lx1*ly1*lz1*1)
1264  $ , h1(lx1,ly1,lz1,1)
1265  $ , h2(lx1,ly1,lz1,1)
1266 
1267  naxhm = naxhm + 1
1268  etime1 = dnekclock()
1269 
1270  nel = nelfld(ifield)
1271  ntot1 = lx1*ly1*lz1*nel
1272 
1273 c if (ifaxis.and.ifsplit) call exitti(
1274 c $'Axisymmetric stress w/PnPn not yet supported.$',istep)
1275 
1276 c icase = 1 --- axsf_fast (no axisymmetry)
1277 c icase = 2 --- stress formulation and supports axisymmetry
1278 c icase = 3 --- 3 separate axhelm calls
1279 
1280  icase = 1 ! Fast mode for stress
1281  if (ifaxis) icase=2 ! Slow for stress, but supports axisymmetry
1282  if (matmod.lt.0) icase=2 ! Elasticity case
1283 c if (matmod.lt.0) icase=3 ! Block-diagonal Axhelm
1284  if (matmod.lt.0) icase=1 ! Elasticity case (faster, 7/28/17,pff)
1285  if (.not.ifstrs) icase=3 ! Block-diagonal Axhelm
1286 
1287  if (icase.eq.1) then
1288 
1289  call axsf_fast(au1,au2,au3,u1,u2,u3,h1,h2,ifield)
1290 
1291  elseif (icase.eq.3) then
1292 
1293  call axhelm(au1,u1,h1,h2,1,1)
1294  call axhelm(au2,u2,h1,h2,1,2)
1295  if (if3d) call axhelm(au3,u3,h1,h2,1,3)
1296 
1297  else ! calculate coupled Aij Uj products
1298 
1299  if ( .not.ifsolv ) call setfast (h1,h2,imesh)
1300 
1301  call stnrate (u1,u2,u3,nel,matmod)
1302  call stress (h1,h2,nel,matmod,ifaxis)
1303  call aijuj (au1,au2,au3,nel,ifaxis)
1304 
1305  if (ifh2 .and. matmod.ge.0) then ! add Helmholtz contributions
1306  call addcol4 (au1,bm1,h2,u1,ntot1)
1307  call addcol4 (au2,bm1,h2,u2,ntot1)
1308  if (ldim.eq.3) call addcol4 (au3,bm1,h2,u3,ntot1)
1309  endif
1310 
1311  endif
1312 
1313  taxhm=taxhm+(dnekclock()-etime1)
1314 
1315  return
1316  end
1317 c-----------------------------------------------------------------------
1318  subroutine stnrate (u1,u2,u3,nel,matmod)
1319 C
1320 C Compute strainrates
1321 C
1322 C CAUTION : Stresses and strainrates share the same scratch commons
1323 C
1324  include 'SIZE'
1325  include 'INPUT'
1326  include 'GEOM'
1327  include 'TSTEP'
1328 
1329  common /ctmp0/ exz(lx1*ly1*lz1*lelt)
1330  $ , eyz(lx1*ly1*lz1*lelt)
1331  common /ctmp1/ exx(lx1*ly1*lz1*lelt)
1332  $ , exy(lx1*ly1*lz1*lelt)
1333  $ , eyy(lx1*ly1*lz1*lelt)
1334  $ , ezz(lx1*ly1*lz1*lelt)
1335 c
1336  dimension u1(lx1,ly1,lz1,1)
1337  $ , u2(lx1,ly1,lz1,1)
1338  $ , u3(lx1,ly1,lz1,1)
1339 C
1340  ntot1 = lx1*ly1*lz1*nel
1341 
1342  CALL rzero3 (exx,eyy,ezz,ntot1)
1343  CALL rzero3 (exy,exz,eyz,ntot1)
1344 
1345  CALL uxyz (u1,exx,exy,exz,nel)
1346  CALL uxyz (u2,exy,eyy,eyz,nel)
1347  IF (ldim.EQ.3) CALL uxyz (u3,exz,eyz,ezz,nel)
1348 
1349  CALL invcol2 (exx,jacm1,ntot1)
1350  CALL invcol2 (exy,jacm1,ntot1)
1351  CALL invcol2 (eyy,jacm1,ntot1)
1352 
1353  IF (ifaxis) CALL axiezz (u2,eyy,ezz,nel)
1354 C
1355  IF (ldim.EQ.3) THEN
1356  CALL invcol2 (exz,jacm1,ntot1)
1357  CALL invcol2 (eyz,jacm1,ntot1)
1358  CALL invcol2 (ezz,jacm1,ntot1)
1359  endif
1360 C
1361  return
1362  end
1363 c-----------------------------------------------------------------------
1364  subroutine stress (h1,h2,nel,matmod,ifaxis)
1365 C
1366 C MATMOD.GE.0 Fluid material models
1367 C MATMOD.LT.0 Solid material models
1368 C
1369 C CAUTION : Stresses and strainrates share the same scratch commons
1370 C
1371  include 'SIZE'
1372  common /ctmp1/ txx(lx1,ly1,lz1,lelt)
1373  $ , txy(lx1,ly1,lz1,lelt)
1374  $ , tyy(lx1,ly1,lz1,lelt)
1375  $ , tzz(lx1,ly1,lz1,lelt)
1376  common /ctmp0/ txz(lx1,ly1,lz1,lelt)
1377  $ , tyz(lx1,ly1,lz1,lelt)
1378  common /scrsf/ t11(lx1,ly1,lz1,lelt)
1379  $ , t22(lx1,ly1,lz1,lelt)
1380  $ , t33(lx1,ly1,lz1,lelt)
1381  $ , hii(lx1,ly1,lz1,lelt)
1382 C
1383  dimension h1(lx1,ly1,lz1,1),h2(lx1,ly1,lz1,1)
1384  LOGICAL IFAXIS
1385 
1386  NTOT1 = lx1*ly1*lz1*nel
1387 
1388  IF (matmod.EQ.0) THEN
1389 
1390 C Newtonian fluids
1391 
1392  const = 2.0
1393  CALL cmult2 (hii,h1,const,ntot1)
1394  CALL col2 (txx,hii,ntot1)
1395  CALL col2 (txy,h1 ,ntot1)
1396  CALL col2 (tyy,hii,ntot1)
1397  IF (ifaxis .OR. ldim.EQ.3) CALL col2 (tzz,hii,ntot1)
1398  IF (ldim.EQ.3) THEN
1399  CALL col2 (txz,h1 ,ntot1)
1400  CALL col2 (tyz,h1 ,ntot1)
1401  endif
1402 C
1403  ELSEIF (matmod.EQ.-1) THEN
1404 C
1405 C Elastic solids
1406 C
1407  const = 2.0
1408  CALL add3s (hii,h1,h2,const,ntot1)
1409  CALL copy (t11,txx,ntot1)
1410  CALL copy (t22,tyy,ntot1)
1411  CALL col3 (txx,hii,t11,ntot1)
1412  CALL addcol3 (txx,h1 ,t22,ntot1)
1413  CALL col3 (tyy,h1 ,t11,ntot1)
1414  CALL addcol3 (tyy,hii,t22,ntot1)
1415  CALL col2 (txy,h2 ,ntot1)
1416  IF (ifaxis .OR. ldim.EQ.3) THEN
1417  CALL copy (t33,tzz,ntot1)
1418  CALL col3 (tzz,h1 ,t11,ntot1)
1419  CALL addcol3 (tzz,h1 ,t22,ntot1)
1420  CALL addcol3 (tzz,hii,t33,ntot1)
1421  CALL addcol3 (txx,h1 ,t33,ntot1)
1422  CALL addcol3 (tyy,h1 ,t33,ntot1)
1423  endif
1424  IF (ldim.EQ.3) THEN
1425  CALL col2 (txz,h2 ,ntot1)
1426  CALL col2 (tyz,h2 ,ntot1)
1427  endif
1428 C
1429  endif
1430 C
1431  return
1432  end
1433 c-----------------------------------------------------------------------
1434  subroutine aijuj (au1,au2,au3,nel,ifaxis)
1435 C
1436  include 'SIZE'
1437  common /ctmp1/ txx(lx1,ly1,lz1,lelt)
1438  $ , txy(lx1,ly1,lz1,lelt)
1439  $ , tyy(lx1,ly1,lz1,lelt)
1440  $ , tzz(lx1,ly1,lz1,lelt)
1441  common /ctmp0/ txz(lx1,ly1,lz1,lelt)
1442  $ , tyz(lx1,ly1,lz1,lelt)
1443 C
1444  dimension au1(lx1,ly1,lz1,1)
1445  $ , au2(lx1,ly1,lz1,1)
1446  $ , au3(lx1,ly1,lz1,1)
1447  LOGICAL IFAXIS
1448 C
1449  CALL ttxyz (au1,txx,txy,txz,nel)
1450  CALL ttxyz (au2,txy,tyy,tyz,nel)
1451  IF (ifaxis) CALL axitzz (au2,tzz,nel)
1452  IF (ldim.EQ.3) CALL ttxyz (au3,txz,tyz,tzz,nel)
1453 C
1454  return
1455  end
1456 c-----------------------------------------------------------------------
1457  subroutine uxyz (u,ex,ey,ez,nel)
1458 C
1459  include 'SIZE'
1460  include 'GEOM'
1461  common /scrsf/ ur(lx1,ly1,lz1,lelt)
1462  $ , us(lx1,ly1,lz1,lelt)
1463  $ , ut(lx1,ly1,lz1,lelt)
1464 c
1465  dimension u(lx1,ly1,lz1,1)
1466  $ , ex(lx1,ly1,lz1,1)
1467  $ , ey(lx1,ly1,lz1,1)
1468  $ , ez(lx1,ly1,lz1,1)
1469 C
1470  ntot1 = lx1*ly1*lz1*nel
1471 C
1472  CALL urst (u,ur,us,ut,nel)
1473 C
1474  CALL addcol3 (ex,rxm1,ur,ntot1)
1475  CALL addcol3 (ex,sxm1,us,ntot1)
1476  CALL addcol3 (ey,rym1,ur,ntot1)
1477  CALL addcol3 (ey,sym1,us,ntot1)
1478 C
1479  IF (ldim.EQ.3) THEN
1480  CALL addcol3 (ez,rzm1,ur,ntot1)
1481  CALL addcol3 (ez,szm1,us,ntot1)
1482  CALL addcol3 (ez,tzm1,ut,ntot1)
1483  CALL addcol3 (ex,txm1,ut,ntot1)
1484  CALL addcol3 (ey,tym1,ut,ntot1)
1485  endif
1486 C
1487  return
1488  end
1489  subroutine urst (u,ur,us,ut,nel)
1490 C
1491  include 'SIZE'
1492  include 'GEOM'
1493  include 'INPUT'
1494 C
1495  dimension u(lx1,ly1,lz1,1)
1496  $ , ur(lx1,ly1,lz1,1)
1497  $ , us(lx1,ly1,lz1,1)
1498  $ , ut(lx1,ly1,lz1,1)
1499 C
1500  DO 100 iel=1,nel
1501  IF (ifaxis) CALL setaxdy ( ifrzer(iel) )
1502  CALL ddrst (u(1,1,1,iel),ur(1,1,1,iel),
1503  $ us(1,1,1,iel),ut(1,1,1,iel))
1504  100 CONTINUE
1505 C
1506  return
1507  end
1508  subroutine ddrst (u,ur,us,ut)
1509 C
1510  include 'SIZE'
1511  include 'DXYZ'
1512 C
1513  dimension u(lx1,ly1,lz1)
1514  $ , ur(lx1,ly1,lz1)
1515  $ , us(lx1,ly1,lz1)
1516  $ , ut(lx1,ly1,lz1)
1517 C
1518  nxy1 = lx1*ly1
1519  nyz1 = ly1*lz1
1520 C
1521  CALL mxm (dxm1,lx1,u,lx1,ur,nyz1)
1522  IF (ldim.EQ.2) THEN
1523  CALL mxm (u,lx1,dytm1,ly1,us,ly1)
1524  ELSE
1525  DO 10 iz=1,lz1
1526  CALL mxm (u(1,1,iz),lx1,dytm1,ly1,us(1,1,iz),ly1)
1527  10 CONTINUE
1528  CALL mxm (u,nxy1,dztm1,lz1,ut,lz1)
1529  endif
1530 C
1531  return
1532  end
1533  subroutine axiezz (u2,eyy,ezz,nel)
1534 C
1535  include 'SIZE'
1536  include 'GEOM'
1537 C
1538  dimension u2(lx1,ly1,lz1,1)
1539  $ , eyy(lx1,ly1,lz1,1)
1540  $ , ezz(lx1,ly1,lz1,1)
1541 C
1542  nxyz1 = lx1*ly1*lz1
1543 C
1544  DO 100 iel=1,nel
1545  IF ( ifrzer(iel) ) THEN
1546  DO 200 ix=1,lx1
1547  ezz(ix, 1,1,iel) = eyy(ix,1,1,iel)
1548  DO 200 iy=2,ly1
1549  ezz(ix,iy,1,iel) = u2(ix,iy,1,iel) / ym1(ix,iy,1,iel)
1550  200 CONTINUE
1551  ELSE
1552  CALL invcol3 (ezz(1,1,1,iel),u2(1,1,1,iel),ym1(1,1,1,iel),
1553  $ nxyz1)
1554  endif
1555  100 CONTINUE
1556 C
1557  return
1558  end
1559 c-----------------------------------------------------------------------
1560  subroutine flush_io
1561  return
1562  end
1563 c-----------------------------------------------------------------------
1564  subroutine fcsum2(xsum,asum,x,e,f)
1565 c
1566 c Compute the weighted sum of X over face f of element e
1567 c
1568 c x is an (NX,NY,NZ) data structure
1569 c f is in the preprocessor notation
1570 c
1571 c xsum is sum (X*area)
1572 c asum is sum (area)
1573 
1574 
1575  include 'SIZE'
1576  include 'GEOM'
1577  include 'TOPOL'
1578  real x(lx1,ly1,lz1,1)
1579  integer e,f,fd
1580 
1581  asum = 0.
1582  xsum = 0.
1583 
1584 c Set up counters ; fd is the dssum notation.
1585  call dsset(lx1,ly1,lz1)
1586  fd = eface1(f)
1587  js1 = skpdat(1,fd)
1588  jf1 = skpdat(2,fd)
1589  jskip1 = skpdat(3,fd)
1590  js2 = skpdat(4,fd)
1591  jf2 = skpdat(5,fd)
1592  jskip2 = skpdat(6,fd)
1593 
1594  i = 0
1595  do j2=js2,jf2,jskip2
1596  do j1=js1,jf1,jskip1
1597  i = i+1
1598  xsum = xsum+area(i,1,f,e)*x(j1,j2,1,e)
1599  asum = asum+area(i,1,f,e)
1600  enddo
1601  enddo
1602 
1603  return
1604  end
1605 c-----------------------------------------------------------------------
1606  function surf_mean(u,ifld,bc_in,ierr)
1607 
1608  include 'SIZE'
1609  include 'TOTAL'
1610 
1611  real u(1)
1612 
1613  integer e,f
1614  character*3 bc_in
1615 
1616  usum = 0
1617  asum = 0
1618 
1619  nface = 2*ldim
1620  do e=1,nelfld(ifld)
1621  do f=1,nface
1622  if (cbc(f,e,ifld).eq.bc_in) then
1623  call fcsum2(usum_f,asum_f,u,e,f)
1624  usum = usum + usum_f
1625  asum = asum + asum_f
1626  endif
1627  enddo
1628  enddo
1629 
1630  usum = glsum(usum,1) ! sum across processors
1631  asum = glsum(asum,1)
1632 
1633  surf_mean = usum
1634  ierr = 1
1635 
1636  if (asum.gt.0) surf_mean = usum/asum
1637  if (asum.gt.0) ierr = 0
1638 
1639  return
1640  end
1641 c-----------------------------------------------------------------------
1642  subroutine fdm_h1a(z,r,d,nel,kt,rr)
1643  include 'SIZE'
1644  include 'TOTAL'
1645 c
1646  common /ctmp0/ w(lx1,ly1,lz1)
1647 c
1648  include 'FDMH1'
1649 c
1650 c Overlapping Schwarz, FDM based
1651 c
1652  real z(lx1,ly1,lz1,1)
1653  real r(lx1,ly1,lz1,1)
1654  real d(lx1,ly1,lz1,1)
1655  real rr(lx1,ly1,lz1,1)
1656 
1657  integer kt(lelt,3)
1658 
1659  integer icalld
1660  save icalld
1661  data icalld /0/
1662 
1663  n1 = lx1
1664  n2 = lx1*lx1
1665  n3 = lx1*lx1*lx1
1666  n = lx1*ly1*lz1*nel
1667 
1668  if (ifbhalf) then
1669  call col3(rr,r,bhalf,n)
1670  else
1671  call copy(rr,r,n)
1672  endif
1673  icalld = icalld+1
1674 c
1675 c
1676  do ie=1,nel
1677  if (if3d) then
1678 c Transfer to wave space:
1679  call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n2)
1680  do iz=1,n1
1681  call mxm(w(1,1,iz),n1,fds(1,kt(ie,2)),n1,z(1,1,iz,ie),n1)
1682  enddo
1683  call mxm(z(1,1,1,ie),n2,fds(1,kt(ie,3)),n1,w,n1)
1684 c
1685 c fdsolve:
1686 c
1687  call col2(w,d(1,1,1,ie),n3)
1688 c
1689 c Transfer to physical space:
1690 c
1691  call mxm(w,n2,fdst(1,kt(ie,3)),n1,z(1,1,1,ie),n1)
1692  do iz=1,n1
1693  call mxm(z(1,1,iz,ie),n1,fdst(1,kt(ie,2)),n1,w(1,1,iz),n1)
1694  enddo
1695  call mxm(fds(1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n2)
1696 c
1697  else
1698 c Transfer to wave space:
1699  call mxm(fdst(1,kt(ie,1)),n1,rr(1,1,1,ie),n1,w,n1)
1700  call mxm(w,n1,fds(1,kt(ie,2)),n1,z(1,1,1,ie),n1)
1701 c
1702 c fdsolve:
1703 c
1704  call col2(z(1,1,1,ie),d(1,1,1,ie),n2)
1705 c
1706 c Transfer to physical space:
1707 c
1708  call mxm(z(1,1,1,ie),n1,fdst(1,kt(ie,2)),n1,w,n1)
1709  call mxm(fds(1,kt(ie,1)),n1,w,n1,z(1,1,1,ie),n1)
1710 c
1711  endif
1712  enddo
1713 
1714  if (ifbhalf) call col2(z,bhalf,n)
1715 
1716  call dssum(z,lx1,ly1,lz1)
1717 
1718  return
1719  end
1720 c-----------------------------------------------------------------------
1721  subroutine set_vert_strs(glo_num,ngv,nx,nel,vertex,ifcenter)
1722 
1723 c Given global array, vertex, pointing to hex vertices, set up
1724 c a new array of global pointers for an nx^ldim set of elements.
1725 
1726  include 'SIZE'
1727  include 'INPUT'
1728 
1729  integer*8 glo_num(1),ngv
1730  integer vertex(1),nx
1731  logical ifcenter
1732 
1733  common /ctmp0/ gnum(lx1*ly1*lz1*lelt)
1734  integer*8 gnum
1735 
1736  call set_vert(gnum,ngv,nx,nel,vertex,ifcenter)
1737 
1738  n=nel*nx*nx
1739  if (if3d) n=n*nx
1740 
1741 c Pack vertex ids component-wise (u,v,w_1; u,v,w_2; etc.)
1742 
1743  k=0
1744  do i=1,n
1745 
1746  do j=1,ldim
1747  glo_num(k+j) = ldim*(gnum(i)-1) + j
1748  enddo
1749  k=k+ldim
1750 
1751  enddo
1752 
1753  return
1754  end
1755 c-----------------------------------------------------------------------
1756  subroutine get_strs_mask(mask,nxc,nzc,nel)
1757  include 'SIZE'
1758  include 'INPUT'
1759  include 'SOLN'
1760 
1761  real mask(ldim,nxc,nxc,nzc,nel)
1762  integer e
1763 
1764 
1765  if (if3d) then
1766  do e=1,nel
1767  do kc=1,nxc
1768  do jc=1,nxc
1769  do ic=1,nxc
1770  k = 1 + ((lx1-1)*(kc-1))/(nxc-1)
1771  j = 1 + ((lx1-1)*(jc-1))/(nxc-1)
1772  i = 1 + ((lx1-1)*(ic-1))/(nxc-1)
1773  mask(1,ic,jc,kc,e) = v1mask(i,j,k,e)
1774  mask(2,ic,jc,kc,e) = v2mask(i,j,k,e)
1775  mask(3,ic,jc,kc,e) = v3mask(i,j,k,e)
1776  enddo
1777  enddo
1778  enddo
1779  enddo
1780  else
1781  do e=1,nel
1782  do jc=1,nxc
1783  do ic=1,nxc
1784  j = 1 + ((lx1-1)*(jc-1))/(nxc-1)
1785  i = 1 + ((lx1-1)*(ic-1))/(nxc-1)
1786  mask(1,ic,jc,1,e) = v1mask(i,j,1,e)
1787  mask(2,ic,jc,1,e) = v2mask(i,j,1,e)
1788  enddo
1789  enddo
1790  enddo
1791  endif
1792 
1793  return
1794  end
1795 c-----------------------------------------------------------------------
1796  subroutine axstrs(a1,a2,a3,p1,p2,p3,h1,h2,matmod,nel)
1797  real a1(1),a2(1),a3(1),p1(1),p2(1),p3(1)
1798 
1799  call axhmsf (a1,a2,a3,p1,p2,p3,h1,h2,matmod)
1800  call rmask (a1,a2,a3,nel)
1801  call opdssum (a1,a2,a3)
1802 
1803  return
1804  end
1805 c-----------------------------------------------------------------------
1806  subroutine axstrs_nds(a1,a2,a3,p1,p2,p3,h1,h2,matmod,nel)
1807  real a1(1),a2(1),a3(1),p1(1),p2(1),p3(1)
1808 
1809  call rmask (p1,p2,p3,nel)
1810  call axhmsf (a1,a2,a3,p1,p2,p3,h1,h2,matmod)
1811  call rmask (a1,a2,a3,nel)
1812 c call opdssum (a1,a2,a3)
1813 
1814  return
1815  end
1816 c-----------------------------------------------------------------------
1817  subroutine get_local_crs_galerkin_strs(a,ncl,nxc,h1,h2,matmod)
1818 
1819 c This routine generates Nelv submatrices of order ncl using
1820 c Galerkin projection
1821 
1822  include 'SIZE'
1823  include 'INPUT'
1824  include 'TSTEP'
1825 
1826 
1827  real a(ldim,ncl,ldim,ncl,1),h1(1),h2(1)
1828 c real a(ncl,ldim,ncl,ldim,1),h1(1),h2(1)
1829 
1830  common /scrcr2/ a1(lx1*ly1*lz1,lelt),w1(lx1*ly1*lz1,lelt)
1831  $ , a2(lx1*ly1*lz1,lelt),w2(lx1*ly1*lz1,lelt)
1832  common /scrcr3/ a3(lx1*ly1*lz1,lelt),w3(lx1*ly1*lz1,lelt)
1833  $ , b(lx1*ly1*lz1,8)
1834 
1835  integer e
1836 
1837  nel = nelfld(ifield)
1838 
1839  do j=1,ncl ! bi- or tri-linear interpolant, ONLY
1840  call gen_crs_basis(b(1,j),j)
1841  enddo
1842 
1843  nxyz = lx1*ly1*lz1
1844  do k = 1,ldim
1845  call opzero(w1,w2,w3)
1846  call opzero(a1,a2,a3)
1847 
1848  do j = 1,ncl
1849 
1850  do e = 1,nel
1851  if (k.eq.1) call copy(w1(1,e),b(1,j),nxyz)
1852  if (k.eq.2) call copy(w2(1,e),b(1,j),nxyz)
1853  if (k.eq.3) call copy(w3(1,e),b(1,j),nxyz)
1854  enddo
1855 
1856  call axstrs_nds(a1,a2,a3,w1,w2,w3,h1,h2,matmod,nel)
1857 
1858  do e = 1,nel
1859  do i = 1,ncl
1860 
1861  a(1,i,k,j,e)=vlsc2(b(1,i),a1(1,e),nxyz) ! bi^T * A^e * bj
1862  a(2,i,k,j,e)=vlsc2(b(1,i),a2(1,e),nxyz) ! bi^T * A^e * bj
1863  if (if3d)
1864  $ a(3,i,k,j,e)=vlsc2(b(1,i),a3(1,e),nxyz) ! bi^T * A^e * bj
1865 
1866  enddo
1867  enddo
1868 
1869  enddo
1870  enddo
1871 
1872  return
1873  end
1874 c-----------------------------------------------------------------------
1875  subroutine crs_strs(u1,u2,u3,v1,v2,v3)
1876 c Given an input vector v, this generates the H1 coarse-grid solution
1877  include 'SIZE'
1878  include 'DOMAIN'
1879  include 'INPUT'
1880  include 'GEOM'
1881  include 'SOLN'
1882  include 'PARALLEL'
1883  include 'CTIMER'
1884  include 'TSTEP'
1885 
1886  real u1(1),u2(1),u3(1),v1(1),v2(1),v3(1)
1887 
1888  common /scrpr3/ uc1(lcr*lelt),uc2(lcr*lelt),uc3(lcr*lelt)
1889  common /scrpr2/ vc1(lcr*lelt),vc2(lcr*lelt),vc3(lcr*lelt)
1890 
1891  integer icalld1
1892  save icalld1
1893  data icalld1 /0/
1894 
1895 
1896  if (icalld1.eq.0) then ! timer info
1897  ncrsl=0
1898  tcrsl=0.0
1899  icalld1=1
1900  endif
1901  ncrsl = ncrsl + 1
1902 
1903  n = nelv*lx1*ly1*lz1
1904  m = nelv*lcr
1905 
1906  call map_f_to_c_h1_bilin(uc1,v1) ! additive Schwarz
1907  call map_f_to_c_h1_bilin(uc2,v2) ! additive Schwarz
1908  if (if3d) call map_f_to_c_h1_bilin(uc3,v3) ! additive Schwarz
1909 
1910  k=0
1911  if (if3d) then
1912  do i=1,m
1913  vc1(k+1)=uc1(i)
1914  vc1(k+2)=uc2(i)
1915  vc1(k+3)=uc3(i)
1916  k=k+3
1917  enddo
1918  else
1919  do i=1,m
1920  vc1(k+1)=uc1(i)
1921  vc1(k+2)=uc2(i)
1922  k=k+2
1923  enddo
1924  endif
1925 
1926  etime1=dnekclock()
1927  call fgslib_crs_solve(xxth_strs,uc1,vc1)
1928  tcrsl=tcrsl+dnekclock()-etime1
1929 
1930  k=0
1931  if (if3d) then
1932  do i=1,m
1933  vc1(i)=uc1(k+1)
1934  vc2(i)=uc1(k+2)
1935  vc3(i)=uc1(k+3)
1936  k=k+3
1937  enddo
1938  else
1939  do i=1,m
1940  vc1(i)=uc1(k+1)
1941  vc2(i)=uc1(k+2)
1942  k=k+2
1943  enddo
1944  endif
1945  call map_c_to_f_h1_bilin(u1,vc1)
1946  call map_c_to_f_h1_bilin(u2,vc2)
1947  if (if3d) call map_c_to_f_h1_bilin(u3,vc3)
1948 
1949  return
1950  end
1951 c-----------------------------------------------------------------------
1952  subroutine set_up_h1_crs_strs(h1,h2,ifld,matmod)
1953 
1954  include 'SIZE'
1955  include 'GEOM'
1956  include 'DOMAIN'
1957  include 'INPUT'
1958  include 'PARALLEL'
1959  include 'TSTEP'
1960  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1961 
1962  common /ivrtx/ vertex((2**ldim)*lelt)
1963  integer vertex
1964 
1965  integer null_space,e
1966 
1967  character*3 cb
1968  common /scrxxti/ ia(ldim*ldim*lcr*lcr*lelv)
1969  $ , ja(ldim*ldim*lcr*lcr*lelv)
1970 
1971  parameter(lcc=2**ldim)
1972  common /scrcr1/ a(ldim*ldim*lcc*lcc*lelt)
1973  real mask(ldim,lcr,lelv)
1974  equivalence (mask,a)
1975 
1976  integer*8 ngv
1977 
1978  t0 = dnekclock()
1979 
1980  nel = nelfld(ifld)
1981 
1982  nxc = 2
1983  nzc = 1
1984  if (if3d) nzc=nxc
1985 
1986  ncr = nxc**ldim
1987  mcr = ldim*ncr ! Dimension of unassembled sub-block
1988  n = mcr*nel
1989 
1990 c Set SEM_to_GLOB index
1991  call get_vertex
1992  call set_vert_strs(se_to_gcrs,ngv,nxc,nel,vertex,.true.)
1993 
1994 
1995 c Set mask for full array
1996  call get_strs_mask (mask,nxc,nzc,nel) ! Set mask
1997  call set_jl_crs_mask (n,mask,se_to_gcrs)
1998 
1999 
2000 c Setup local SEM-based Neumann operators (for now, just full...)
2001  call get_local_crs_galerkin_strs(a,ncr,nxc,h1,h2,matmod)
2002  call set_mat_ij(ia,ja,mcr,nel)
2003  nnz=mcr*mcr*nel
2004 
2005  null_space=0
2006 
2007  t1 = dnekclock()-t0
2008  if (nio.eq.0)
2009  $ write(6,*) 'start:: setup h1 coarse grid ',t1, ' sec'
2010  $ ,nnz,mcr,ncr,nel
2011 
2012  do i=1,nnz
2013  write(44,44) ia(i),ja(i),a(i)
2014  enddo
2015  44 format(2i9,1pe22.13)
2016 c stop
2017 
2018  imode = param(40)
2019  call fgslib_crs_setup(xxth_strs,imode,nekcomm,mp,n,se_to_gcrs,
2020  $ nnz,ia,ja,a,null_space)
2021 
2022  t0 = dnekclock()-t0
2023  if (nio.eq.0) then
2024  write(6,*) 'done :: setup h1 coarse grid ',t0, ' sec',xxth_strs
2025  write(6,*) ' '
2026  endif
2027 
2028  return
2029  end
2030 c-----------------------------------------------------------------------
2031  subroutine axsf_e_3d(au,av,aw,u,v,w,h1,h2,ur,e)
2032 c
2033 c du_i
2034 c Compute the gradient tensor G_ij := ---- , for element e
2035 c du_j
2036 c
2037  include 'SIZE'
2038  include 'TOTAL'
2039 
2040  real au(1),av(1),aw(1),u(1),v(1),w(1),h1(1),h2(1)
2041 
2042  parameter(l=lx1*ly1*lz1)
2043  real ur(l,ldim,ldim)
2044 
2045  integer e,p
2046 
2047  p = lx1-1 ! Polynomial degree
2048  n = lx1*ly1*lz1
2049 
2050 
2051  call local_grad3(ur(1,1,1),ur(1,2,1),ur(1,3,1),u,p,1,dxm1,dxtm1)
2052  call local_grad3(ur(1,1,2),ur(1,2,2),ur(1,3,2),v,p,1,dxm1,dxtm1)
2053  call local_grad3(ur(1,1,3),ur(1,2,3),ur(1,3,3),w,p,1,dxm1,dxtm1)
2054 
2055  do i=1,n
2056 
2057  u1 = ur(i,1,1)*rxm1(i,1,1,e) + ur(i,2,1)*sxm1(i,1,1,e)
2058  $ + ur(i,3,1)*txm1(i,1,1,e)
2059  u2 = ur(i,1,1)*rym1(i,1,1,e) + ur(i,2,1)*sym1(i,1,1,e)
2060  $ + ur(i,3,1)*tym1(i,1,1,e)
2061  u3 = ur(i,1,1)*rzm1(i,1,1,e) + ur(i,2,1)*szm1(i,1,1,e)
2062  $ + ur(i,3,1)*tzm1(i,1,1,e)
2063 
2064  v1 = ur(i,1,2)*rxm1(i,1,1,e) + ur(i,2,2)*sxm1(i,1,1,e)
2065  $ + ur(i,3,2)*txm1(i,1,1,e)
2066  v2 = ur(i,1,2)*rym1(i,1,1,e) + ur(i,2,2)*sym1(i,1,1,e)
2067  $ + ur(i,3,2)*tym1(i,1,1,e)
2068  v3 = ur(i,1,2)*rzm1(i,1,1,e) + ur(i,2,2)*szm1(i,1,1,e)
2069  $ + ur(i,3,2)*tzm1(i,1,1,e)
2070 
2071  w1 = ur(i,1,3)*rxm1(i,1,1,e) + ur(i,2,3)*sxm1(i,1,1,e)
2072  $ + ur(i,3,3)*txm1(i,1,1,e)
2073  w2 = ur(i,1,3)*rym1(i,1,1,e) + ur(i,2,3)*sym1(i,1,1,e)
2074  $ + ur(i,3,3)*tym1(i,1,1,e)
2075  w3 = ur(i,1,3)*rzm1(i,1,1,e) + ur(i,2,3)*szm1(i,1,1,e)
2076  $ + ur(i,3,3)*tzm1(i,1,1,e)
2077 
2078  dj = h1(i)*w3m1(i,1,1)*jacmi(i,e)
2079  s11 = dj*(u1 + u1)! S_ij
2080  s12 = dj*(u2 + v1)
2081  s13 = dj*(u3 + w1)
2082  s21 = dj*(v1 + u2)
2083  s22 = dj*(v2 + v2)
2084  s23 = dj*(v3 + w2)
2085  s31 = dj*(w1 + u3)
2086  s32 = dj*(w2 + v3)
2087  s33 = dj*(w3 + w3)
2088 
2089 c Sum_j : (r_p/x_j) h1 J S_ij
2090 
2091  ur(i,1,1)=rxm1(i,1,1,e)*s11+rym1(i,1,1,e)*s12+rzm1(i,1,1,e)*s13
2092  ur(i,2,1)=sxm1(i,1,1,e)*s11+sym1(i,1,1,e)*s12+szm1(i,1,1,e)*s13
2093  ur(i,3,1)=txm1(i,1,1,e)*s11+tym1(i,1,1,e)*s12+tzm1(i,1,1,e)*s13
2094 
2095  ur(i,1,2)=rxm1(i,1,1,e)*s21+rym1(i,1,1,e)*s22+rzm1(i,1,1,e)*s23
2096  ur(i,2,2)=sxm1(i,1,1,e)*s21+sym1(i,1,1,e)*s22+szm1(i,1,1,e)*s23
2097  ur(i,3,2)=txm1(i,1,1,e)*s21+tym1(i,1,1,e)*s22+tzm1(i,1,1,e)*s23
2098 
2099  ur(i,1,3)=rxm1(i,1,1,e)*s31+rym1(i,1,1,e)*s32+rzm1(i,1,1,e)*s33
2100  ur(i,2,3)=sxm1(i,1,1,e)*s31+sym1(i,1,1,e)*s32+szm1(i,1,1,e)*s33
2101  ur(i,3,3)=txm1(i,1,1,e)*s31+tym1(i,1,1,e)*s32+tzm1(i,1,1,e)*s33
2102 
2103  enddo
2104  call local_grad3_t
2105  $ (au,ur(1,1,1),ur(1,2,1),ur(1,3,1),p,1,dxm1,dxtm1,av)
2106  call local_grad3_t
2107  $ (av,ur(1,1,2),ur(1,2,2),ur(1,3,2),p,1,dxm1,dxtm1,ur)
2108  call local_grad3_t
2109  $ (aw,ur(1,1,3),ur(1,2,3),ur(1,3,3),p,1,dxm1,dxtm1,ur)
2110 
2111  do i=1,n
2112  au(i)=au(i) + h2(i)*bm1(i,1,1,e)*u(i)
2113  av(i)=av(i) + h2(i)*bm1(i,1,1,e)*v(i)
2114  aw(i)=aw(i) + h2(i)*bm1(i,1,1,e)*w(i)
2115  enddo
2116 
2117  return
2118  end
2119 c-----------------------------------------------------------------------
2120  subroutine axsf_e_2d(au,av,u,v,h1,h2,ur,e)
2121 c
2122 c du_i
2123 c Compute the gradient tensor G_ij := ---- , for element e
2124 c du_j
2125 c
2126  include 'SIZE'
2127  include 'TOTAL'
2128 
2129  real au(1),av(1),u(1),v(1),h1(1),h2(1)
2130 
2131  parameter(l=lx1*ly1*lz1)
2132  real ur(l,ldim,ldim)
2133 
2134  integer e,p
2135 
2136  p = lx1-1 ! Polynomial degree
2137  n = lx1*ly1*lz1
2138 
2139 
2140  call local_grad2(ur(1,1,1),ur(1,2,1),u,p,1,dxm1,dxtm1)
2141  call local_grad2(ur(1,1,2),ur(1,2,2),v,p,1,dxm1,dxtm1)
2142 
2143  do i=1,n
2144  dj = h1(i)*w3m1(i,1,1)*jacmi(i,e)
2145 
2146  u1 = ur(i,1,1)*rxm1(i,1,1,e) + ur(i,2,1)*sxm1(i,1,1,e) !ux
2147  u2 = ur(i,1,1)*rym1(i,1,1,e) + ur(i,2,1)*sym1(i,1,1,e) !uy
2148  v1 = ur(i,1,2)*rxm1(i,1,1,e) + ur(i,2,2)*sxm1(i,1,1,e) !vx
2149  v2 = ur(i,1,2)*rym1(i,1,1,e) + ur(i,2,2)*sym1(i,1,1,e) !vy
2150 
2151  s11 = dj*( u1 + u1 ) ! h1*rho*S_ij
2152  s12 = dj*( u2 + v1 )
2153  s21 = dj*( v1 + u2 )
2154  s22 = dj*( v2 + v2 )
2155 
2156 c Sum_j : (r_k/x_j) h1 J S_ij
2157 
2158  ur(i,1,1)=rxm1(i,1,1,e)*s11+rym1(i,1,1,e)*s12 ! i=1,k=1
2159  ur(i,2,1)=sxm1(i,1,1,e)*s11+sym1(i,1,1,e)*s12 ! i=1,k=2
2160 
2161  ur(i,1,2)=rxm1(i,1,1,e)*s21+rym1(i,1,1,e)*s22 ! i=2,k=1
2162  ur(i,2,2)=sxm1(i,1,1,e)*s21+sym1(i,1,1,e)*s22 ! i=2,k=2
2163 
2164  enddo
2165 
2166  call local_grad2_t(au,ur(1,1,1),ur(1,2,1),p,1,dxm1,dxtm1,av)
2167  call local_grad2_t(av,ur(1,1,2),ur(1,2,2),p,1,dxm1,dxtm1,ur)
2168 
2169  do i=1,n
2170  au(i)=au(i) + h2(i)*bm1(i,1,1,e)*u(i)
2171  av(i)=av(i) + h2(i)*bm1(i,1,1,e)*v(i)
2172  enddo
2173 
2174  return
2175  end
2176 c-----------------------------------------------------------------------
2177  subroutine axsf_fast(au,av,aw,u,v,w,h1,h2,ifld)
2178  include 'SIZE'
2179  include 'TOTAL'
2180 
2181  parameter(l=lx1*ly1*lz1)
2182  real au(l,1),av(l,1),aw(l,1),u(l,1),v(l,1),w(l,1),h1(l,1),h2(l,1)
2183 
2184  common /btmp0/ ur(l,ldim,ldim)
2185 
2186  integer e
2187 
2188  nel = nelfld(ifld)
2189 
2190  if (if3d) then
2191  do e=1,nel
2192  call axsf_e_3d(au(1,e),av(1,e),aw(1,e),u(1,e),v(1,e),w(1,e)
2193  $ ,h1(1,e),h2(1,e),ur,e)
2194  enddo
2195  else
2196  do e=1,nel
2197  call axsf_e_2d(au(1,e),av(1,e),u(1,e),v(1,e)
2198  $ ,h1(1,e),h2(1,e),ur,e)
2199  enddo
2200  endif
2201 
2202  return
2203  end
2204 c-----------------------------------------------------------------------
2205  subroutine ttxyz (ff,tx,ty,tz,nel)
2206 C
2207  include 'SIZE'
2208  include 'DXYZ'
2209  include 'GEOM'
2210  include 'INPUT'
2211  include 'TSTEP'
2212  include 'WZ'
2213 
2214  dimension tx(lx1,ly1,lz1,1)
2215  $ , ty(lx1,ly1,lz1,1)
2216  $ , tz(lx1,ly1,lz1,1)
2217  $ , ff(lx1*ly1*lz1,1)
2218 
2219  common /scrsf/ fr(lx1*ly1*lz1,lelt)
2220  $ , fs(lx1*ly1*lz1,lelt)
2221  $ , ft(lx1*ly1*lz1,lelt)
2222  real wa(lx1,ly1,lz1,lelt)
2223  equivalence(wa,ft)
2224  real ys(lx1)
2225 
2226  nxyz1 = lx1*ly1*lz1
2227  ntot1 = nxyz1*nel
2228 
2229  CALL col3 (fr,rxm1,tx,ntot1)
2230  CALL addcol3 (fr,rym1,ty,ntot1)
2231  CALL col3 (fs,sxm1,tx,ntot1)
2232  CALL addcol3 (fs,sym1,ty,ntot1)
2233 
2234  IF (ldim.EQ.3) THEN
2235  CALL addcol3 (fr,rzm1,tz,ntot1)
2236  CALL addcol3 (fs,szm1,tz,ntot1)
2237  CALL col3 (ft,txm1,tx,ntot1)
2238  CALL addcol3 (ft,tym1,ty,ntot1)
2239  CALL addcol3 (ft,tzm1,tz,ntot1)
2240  endif
2241 C
2242  IF (ifaxis) THEN
2243  DO 100 iel=1,nel
2244  IF ( ifrzer(iel) ) THEN
2245  CALL mxm (ym1(1,1,1,iel),lx1,datm1,ly1,ys,1)
2246  DO 120 ix=1,lx1
2247  iy = 1
2248  wa(ix,iy,1,iel)=ys(ix)*w2am1(ix,iy)
2249  DO 120 iy=2,ly1
2250  dnr = 1.0 + zam1(iy)
2251  wa(ix,iy,1,iel)=ym1(ix,iy,1,iel)*w2am1(ix,iy)/dnr
2252  120 CONTINUE
2253  ELSE
2254  CALL col3 (wa(1,1,1,iel),ym1(1,1,1,iel),w2cm1,nxyz1)
2255  endif
2256  100 CONTINUE
2257  CALL col2 (fr,wa,ntot1)
2258  CALL col2 (fs,wa,ntot1)
2259  else
2260  do 180 iel=1,nel
2261  call col2(fr(1,iel),w3m1,nxyz1)
2262  call col2(fs(1,iel),w3m1,nxyz1)
2263  call col2(ft(1,iel),w3m1,nxyz1)
2264  180 continue
2265  endif
2266 
2267 
2268  DO 200 iel=1,nel
2269  IF (ifaxis) CALL setaxdy ( ifrzer(iel) )
2270  CALL ttrst (ff(1,iel),fr(1,iel),fs(1,iel),
2271  $ ft(1,iel),fr(1,iel)) ! FR work array
2272  200 CONTINUE
2273 C
2274  return
2275  end
2276 c-----------------------------------------------------------------------
2277  subroutine ttrst (ff,fr,fs,ft,ta)
2278 
2279  include 'SIZE'
2280  include 'DXYZ'
2281  include 'TSTEP'
2282 
2283  dimension ff(lx1,ly1,lz1)
2284  $ , fr(lx1,ly1,lz1)
2285  $ , fs(lx1,ly1,lz1)
2286  $ , ft(lx1,ly1,lz1)
2287  $ , ta(lx1,ly1,lz1)
2288 
2289  nxy1 = lx1*ly1
2290  nyz1 = ly1*lz1
2291  nxyz1 = nxy1*lz1
2292 
2293  CALL mxm (dxtm1,lx1,fr,lx1,ff,nyz1)
2294  IF (ldim.EQ.2) THEN
2295  CALL mxm (fs,lx1,dym1,ly1,ta,ly1)
2296  CALL add2 (ff,ta,nxyz1)
2297  ELSE
2298  DO 10 iz=1,lz1
2299  CALL mxm (fs(1,1,iz),lx1,dym1,ly1,ta(1,1,iz),ly1)
2300  10 CONTINUE
2301  CALL add2 (ff,ta,nxyz1)
2302  CALL mxm (ft,nxy1,dzm1,lz1,ta,lz1)
2303  CALL add2 (ff,ta,nxyz1)
2304  endif
2305 
2306  return
2307  end
2308 c-----------------------------------------------------------------------
2309  subroutine axitzz (vfy,tzz,nel)
2310 C
2311  include 'SIZE'
2312  include 'DXYZ'
2313  include 'GEOM'
2314  include 'WZ'
2315  common /ctmp0/ phi(lx1,ly1)
2316 C
2317  dimension vfy(lx1,ly1,lz1,1)
2318  $ , tzz(lx1,ly1,lz1,1)
2319 C
2320  nxyz1 = lx1*ly1*lz1
2321 C
2322  DO 100 iel=1,nel
2323  CALL setaxw1 ( ifrzer(iel) )
2324  CALL col4 (phi,tzz(1,1,1,iel),jacm1(1,1,1,iel),w3m1,nxyz1)
2325  IF ( ifrzer(iel) ) THEN
2326  DO 120 ix=1,lx1
2327  DO 120 iy=2,ly1
2328  dnr = phi(ix,iy)/( 1.0 + zam1(iy) )
2329  dds = wxm1(ix) * wam1(1) * datm1(iy,1) *
2330  $ jacm1(ix,1,1,iel) * tzz(ix,1,1,iel)
2331  vfy(ix,iy,1,iel)=vfy(ix,iy,1,iel) + dnr + dds
2332  120 CONTINUE
2333  ELSE
2334  CALL add2 (vfy(1,1,1,iel),phi,nxyz1)
2335  endif
2336  100 CONTINUE
2337 C
2338  return
2339  end
2340 c-----------------------------------------------------------------------
2341  subroutine setaxdy (ifaxdy)
2342 C
2343  include 'SIZE'
2344  include 'DXYZ'
2345 C
2346  LOGICAL IFAXDY
2347 C
2348  IF (ifaxdy) THEN
2349  CALL copy (dym1 ,dam1 ,ly1*ly1)
2350  CALL copy (dytm1,datm1,ly1*ly1)
2351  ELSE
2352  CALL copy (dym1 ,dcm1 ,ly1*ly1)
2353  CALL copy (dytm1,dctm1,ly1*ly1)
2354  endif
2355 C
2356  return
2357  end
2358 c-----------------------------------------------------------------------
2359  function opnorm2w(v1,v2,v3,w)
2360  include 'SIZE'
2361  include 'TOTAL'
2362 c
2363  real v1(1) , v2(1), v3(1), w(1)
2364 
2365  n=lx1*ly1*lz1*nelv
2366  s=vlsc3(v1,w,v1,n)+vlsc3(v2,w,v2,n)
2367  if(if3d) s=s+vlsc3(v3,w,v3,n)
2368  s=glsum(s,1)
2369 
2370  if (s.gt.0) s=sqrt(s/volvm1)
2371  opnorm2w = s
2372 
2373  return
2374  end
2375 c-----------------------------------------------------------------------
2376  subroutine strs_project_a(b1,b2,b3,h1,h2,wt,ifld,ierr,matmod)
2377 
2378 c Assumes if uservp is true and thus reorthogonalizes every step
2379 
2380  include 'SIZE'
2381  include 'TOTAL'
2382  include 'ORTHOSTRS' ! Storage of approximation space
2383  include 'CTIMER'
2384 
2385  real b1(1),b2(1),b3(1),h1(1),h2(1),wt(1)
2386  common /ctmp1/ w(lx1*ly1*lz1*lelt,ldim)
2387  real l2a,l2b
2388 
2389  kmax = napproxstrs(1)
2390  k = napproxstrs(2)
2391  n = lx1*ly1*lz1*nelv
2392  m = n*ldim
2393 
2394  if (k.eq.0.or.kmax.eq.0) return
2395 
2396  etime0 = dnekclock()
2397 
2398  l2b=opnorm2w(b1,b2,b3,binvm1)
2399 
2400 c Reorthogonalize basis
2401 
2402  dh1max=difmax(bstrs(1 ),h1,n)
2403  call col3(bstrs,h2,bm1,n)
2404  dh2max=difmax(bstrs(1+n),bstrs,n)
2405 
2406  if (dh1max.gt.0.or.dh2max.gt.0) then
2407  call strs_ortho_all(xstrs(1+m),bstrs(1+m),n,k,h1,h2,wt,ifld,w
2408  $ ,ierr,matmod)
2409  else
2410  call strs_ortho_one(xstrs(1+m),bstrs(1+m),n,k,h1,h2,wt,ifld,w
2411  $ ,ierr,matmod)
2412  endif
2413 
2414  napproxstrs(2) = k
2415 
2416  call opcopy(bstrs(1),bstrs(1+n),bstrs(1+2*n),b1,b2,b3)
2417  call opzero(xstrs(1),xstrs(1+n),xstrs(1+2*n))
2418 
2419  do i=1,k
2420  i1 = 1 + 0*n + (i-1)*m + m
2421  i2 = 1 + 1*n + (i-1)*m + m
2422  i3 = 1 + 2*n + (i-1)*m + m
2423  alpha=op_glsc2_wt(bstrs(1),bstrs(1+n),bstrs(1+2*n)
2424  $ ,xstrs(i1),xstrs(i2),xstrs(i3),wt)
2425 
2426  alphm=-alpha
2427  call opadds(bstrs(1),bstrs(1+n),bstrs(1+2*n)
2428  $ ,bstrs(i1),bstrs(i2),bstrs(i3),alphm,n,2)
2429 
2430  call opadds(xstrs(1),xstrs(1+n),xstrs(1+2*n)
2431  $ ,xstrs(i1),xstrs(i2),xstrs(i3),alpha,n,2)
2432 
2433 
2434  enddo
2435 
2436  call opcopy(b1,b2,b3,bstrs(1),bstrs(1+n),bstrs(1+2*n))
2437  l2a=opnorm2w(b1,b2,b3,binvm1)
2438 
2439  call copy(bstrs(1 ),h1,n) ! Save h1 and h2 for comparison
2440  call col3(bstrs(1+n),bm1,h2,n)
2441 
2442  ratio=l2b/l2a
2443  if (nio.eq.0) write(6,1) istep,' Project ' // 'HELM3 ',
2444  $ l2a,l2b,ratio,k,kmax
2445  1 format(i11,a,6x,1p3e13.4,i4,i4)
2446 
2447 c if (ierr.ne.0) call exitti(' h3proj quit$',ierr)
2448 
2449  tproj = tproj + dnekclock() - etime0
2450 
2451  return
2452  end
2453 c-----------------------------------------------------------------------
2454  subroutine strs_project_b(x1,x2,x3,h1,h2,wt,ifld,ierr)
2455 
2456 c Reconstruct solution; don't bother to orthonomalize bases
2457 
2458  include 'SIZE'
2459  include 'TOTAL'
2460  include 'ORTHOSTRS' ! Storage of approximation space
2461  include 'CTIMER'
2462 
2463  real x1(1),x2(1),x3(1),h1(1),h2(1),wt(1)
2464  common /cptst/ xs(lx1*ly1*lz1*lelt*ldim)
2465 
2466  kmax = napproxstrs(1)
2467  k = napproxstrs(2)
2468  n = lx1*ly1*lz1*nelv
2469  m = n*ldim
2470 
2471  if (kmax.eq.0) return
2472 
2473  etime0 = dnekclock()
2474 
2475  if (k.eq.0) then ! _
2476  call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx+x
2477  k=1
2478  k1 = 1 + 0*n + (k-1)*m + m
2479  k2 = 1 + 1*n + (k-1)*m + m
2480  k3 = 1 + 2*n + (k-1)*m + m
2481  call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! x1=x^n
2482  elseif (k.eq.kmax) then ! _
2483  call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx+x
2484  k=1
2485  k1 = 1 + 0*n + (k-1)*m + m
2486  k2 = 1 + 1*n + (k-1)*m + m
2487  k3 = 1 + 2*n + (k-1)*m + m
2488  call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! x1=x^n
2489 c k=2
2490 c k1 = 1 + 0*n + (k-1)*m + m
2491 c k2 = 1 + 1*n + (k-1)*m + m
2492 c k3 = 1 + 2*n + (k-1)*m + m
2493 c call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),xs(1),xs(1+n),xs(1+2*n))
2494  else
2495  k=k+1
2496  k1 = 1 + 0*n + (k-1)*m + m
2497  k2 = 1 + 1*n + (k-1)*m + m
2498  k3 = 1 + 2*n + (k-1)*m + m
2499  call opcopy(xstrs(k1),xstrs(k2),xstrs(k3),x1,x2,x3) ! xk=dx _
2500  call opadd2(x1,x2,x3,xstrs(1),xstrs(1+n),xstrs(1+2*n)) ! x=dx + x
2501  endif
2502 
2503 c if (k.eq.kmax) call opcopy(xs(1),xs(1+n),xs(1+2*n),x1,x2,x3) ! presave
2504 
2505  napproxstrs(2)=k
2506 
2507  tproj = tproj + dnekclock() - etime0
2508 
2509  return
2510  end
2511 c-----------------------------------------------------------------------
2512  subroutine strs_orthok(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2513 
2514 c Orthonormalize the kth element of X against x_j, j < k.
2515 
2516  include 'SIZE'
2517  include 'TOTAL'
2518 
2519  real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2520  real al(mxprev),bt(mxprev)
2521 
2522  m = n*ldim ! vector length
2523  nel = nelfld(ifld)
2524 
2525  call axhmsf (b(1,1,k),b(1,2,k),b(1,3,k)
2526  $ ,x(1,1,k),x(1,2,k),x(1,3,k),h1,h2,matmod)
2527  call rmask (b(1,1,k),b(1,2,k),b(1,3,k),nel)
2528  call opdssum (b(1,1,k),b(1,2,k),b(1,3,k))
2529 
2530  xax0 = op_glsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2531  $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2532  ierr = 3
2533 
2534  if (xax0.le.0.and.nio.eq.0)write(6,*)istep,ierr,k,xax0,'Proj3 ERR'
2535  if (xax0.le.0) return
2536 
2537  s = 0.
2538  do j=1,k-1! Modifed Gram-Schmidt
2539 
2540  betaj = ( op_vlsc2_wt(b(1,1,j),b(1,2,j),b(1,3,j)
2541  $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2542  $ + op_vlsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2543  $ ,x(1,1,j),x(1,2,j),x(1,3,j),wt))/2.
2544  betam = -glsum(betaj,1)
2545  call add2s2 (x(1,1,k),x(1,1,j),betam,m) ! Full-vector-subtract
2546  call add2s2 (b(1,1,k),b(1,1,j),betam,m) !
2547 
2548  s = s + betam**2
2549 
2550  enddo
2551 
2552  xax1 = xax0-s
2553  xax2 = op_glsc2_wt(b(1,1,k),b(1,2,k),b(1,3,k)
2554  $ ,x(1,1,k),x(1,2,k),x(1,3,k),wt)
2555  scale = xax2
2556 
2557  eps = 1.e-8
2558  ierr = 0
2559  if (scale/xax0.lt.eps) ierr=1
2560 
2561 c if(nio.eq.0) write(6,*) istep,ierr,k,scale,xax0,' SCALE'
2562 c if(nio.eq.0.and.(istep.lt.10.or.mod(istep,100).eq.0.or.ierr.gt.0))
2563 c $write(6,3) istep,k,ierr,xax1/xax0,xax2/xax0
2564 c 3 format(i9,2i3,1p2e12.4,' scale ortho')
2565 
2566  if (scale.gt.0) then
2567  scale = 1./sqrt(scale)
2568  call cmult(x(1,1,k),scale,m)
2569  call cmult(b(1,1,k),scale,m)
2570  else
2571  ierr=2
2572  endif
2573 
2574  return
2575  end
2576 c-----------------------------------------------------------------------
2577  subroutine strs_ortho_one(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2578 
2579  include 'SIZE'
2580  include 'TOTAL'
2581 
2582  real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2583 
2584  m = n*ldim
2585 
2586  js=k
2587  do j=k,k
2588  if (js.lt.j) call copy(x(1,1,js),x(1,1,j),m)
2589  call strs_orthok(x,b,n,js,h1,h2,wt,ifld,w,ierr,matmod)
2590  if (ierr.eq.0) js=js+1
2591  enddo
2592  k = js-1
2593 
2594  return
2595  end
2596 c-----------------------------------------------------------------------
2597  subroutine strs_ortho_all(x,b,n,k,h1,h2,wt,ifld,w,ierr,matmod)
2598 
2599  include 'SIZE'
2600  include 'TOTAL'
2601 
2602  real x(n,ldim,k),b(n,ldim,k),h1(n),h2(n),wt(n),w(n,ldim)
2603 
2604  m = n*ldim
2605 
2606  js=1
2607  do j=1,k
2608  if (js.lt.j) call copy(x(1,1,js),x(1,1,j),m)
2609  call strs_orthok(x,b,n,js,h1,h2,wt,ifld,w,ierr,matmod)
2610  if (ierr.eq.0) js=js+1
2611  enddo
2612  k = js-1
2613 
2614  return
2615  end
2616 c-----------------------------------------------------------------------
2617  subroutine setprop
2618 C------------------------------------------------------------------------
2619 C
2620 C Set variable property arrays
2621 C
2622 C------------------------------------------------------------------------
2623  include 'SIZE'
2624  include 'TOTAL'
2625  include 'CTIMER'
2626 C
2627 C Caution: 2nd and 3rd strainrate invariants residing in scratch
2628 C common /SCREV/ are used in STNRINV and NEKASGN
2629 C
2630  common /screv/ sii(lx1,ly1,lz1,lelt),siii(lx1,ly1,lz1,lelt)
2631 
2632  if (nio.eq.0.and.loglevel.gt.2)
2633  $ write(6,*) 'setprop'
2634 
2635 #ifdef TIMER
2636  if (icalld.eq.0) tspro=0.0
2637  icalld=icalld+1
2638  nspro=icalld
2639  etime1=dnekclock()
2640 #endif
2641 
2642  nxyz1 = lx1*ly1*lz1
2643  mfield=2
2644  IF (ifflow) mfield=1
2645  nfldt = nfield
2646  if (ifmhd) nfldt = nfield+1
2647 
2648  ifld = ifield
2649 
2650  do ifield=mfield,nfldt
2651  if (idpss(ifield-1).eq.-1) goto 100
2652  call vprops
2653  100 continue
2654  enddo
2655 
2656  ifield = ifld
2657 
2658 #ifdef TIMER
2659  tspro=tspro+(dnekclock()-etime1)
2660 #endif
2661 
2662  return
2663  end
2664 c-----------------------------------------------------------------------
setdrs
subroutine setdrs(DRM1, DRTM1, DSM1, DSTM1, IFC)
Definition: bdry.f:1692
faceis
subroutine faceis(CB, S, IEL, IFACE, NX, NY, NZ)
Definition: bdry.f:919
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
nekasgn
subroutine nekasgn(ix, iy, iz, e)
Definition: bdry.f:1121
rzero3
subroutine rzero3(A, B, C, N)
Definition: bdry.f:1821
facec2
subroutine facec2(A1, A2, B1, B2, IFC)
Definition: bdry.f:1782
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
setfast
subroutine setfast(helm1, helm2, imesh)
Definition: hmholtz.f:263
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
set_fdm_prec_h1b
subroutine set_fdm_prec_h1b(d, h1, h2, nel)
Definition: hmholtz.f:1274
set_fdm_prec_h1a
subroutine set_fdm_prec_h1a
Definition: hmholtz.f:1360
setprec
subroutine setprec(dpcm1, helm1, helm2, imsh, isd)
Definition: hmholtz.f:381
opzero
subroutine opzero(ux, uy, uz)
Definition: induct.f:406
compute_cfl
subroutine compute_cfl(cfl, u, v, w, dt)
Definition: induct.f:918
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
glmin
function glmin(a, n)
Definition: math.f:973
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
difmax
real function difmax(a, b, n)
Definition: math.f:1618
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
glmin_ms
function glmin_ms(a, n)
Definition: math.f:1633
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
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
add4
subroutine add4(a, b, c, d, n)
Definition: math.f:728
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
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
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
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
makeuf
subroutine makeuf
Definition: navier1.f:1464
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
local_grad3_t
subroutine local_grad3_t(u, ur, us, ut, N, e, D, Dt, w)
Definition: navier1.f:4669
opcopy
subroutine opcopy(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2474
opcolv
subroutine opcolv(a1, a2, a3, c)
Definition: navier1.f:2418
normvc
subroutine normvc(h1, semi, l2, linf, x1, x2, x3)
Definition: navier1.f:2083
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opsub3
subroutine opsub3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
Definition: navier1.f:2370
local_grad2_t
subroutine local_grad2_t(u, ur, us, N, e, D, Dt, w)
Definition: navier1.f:4694
normsc
subroutine normsc(h1, semi, l2, linf, x, imesh)
Definition: navier1.f:2026
vlsc3
function vlsc3(X, Y, B, N)
Definition: navier4.f:324
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
map_c_to_f_h1_bilin
subroutine map_c_to_f_h1_bilin(uf, uc)
Definition: navier8.f:1407
get_vertex
subroutine get_vertex
Definition: navier8.f:1629
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4
map_f_to_c_h1_bilin
subroutine map_f_to_c_h1_bilin(uc, uf)
Definition: navier8.f:1456
gen_crs_basis
subroutine gen_crs_basis(b, j)
Definition: navier8.f:1540
set_mat_ij
subroutine set_mat_ij(ia, ja, n, ne)
Definition: navier8.f:247
set_jl_crs_mask
subroutine set_jl_crs_mask(n, mask, se_to_gcrs)
Definition: navier8.f:238
faccl2
subroutine faccl2(a, b, iface1)
Definition: subs1.f:909
strs_project_b
subroutine strs_project_b(x1, x2, x3, h1, h2, wt, ifld, ierr)
Definition: subs1.f:2455
ttrst
subroutine ttrst(ff, fr, fs, ft, ta)
Definition: subs1.f:2278
setdt
subroutine setdt
Definition: subs1.f:179
ddrst
subroutine ddrst(u, ur, us, ut)
Definition: subs1.f:1509
chktcgs
subroutine chktcgs(r1, r2, r3, rmask1, rmask2, rmask3, rmult, binv, vol, tol, nel)
Definition: subs1.f:1160
axitzz
subroutine axitzz(vfy, tzz, nel)
Definition: subs1.f:2310
faddcl3
subroutine faddcl3(a, b, c, iface1)
Definition: subs1.f:978
flush_io
subroutine flush_io
Definition: subs1.f:1561
strs_project_a
subroutine strs_project_a(b1, b2, b3, h1, h2, wt, ifld, ierr, matmod)
Definition: subs1.f:2377
axsf_e_3d
subroutine axsf_e_3d(au, av, aw, u, v, w, h1, h2, ur, e)
Definition: subs1.f:2032
nekuvp
subroutine nekuvp(iel)
Definition: subs1.f:1053
fcsum2
subroutine fcsum2(xsum, asum, x, e, f)
Definition: subs1.f:1565
axhmsf
subroutine axhmsf(au1, au2, au3, u1, u2, u3, h1, h2, matmod)
Definition: subs1.f:1241
setdtfs
subroutine setdtfs(dtfs)
Definition: subs1.f:680
set_up_h1_crs_strs
subroutine set_up_h1_crs_strs(h1, h2, ifld, matmod)
Definition: subs1.f:1953
cvgnlps
subroutine cvgnlps(ifconv)
Definition: subs1.f:324
cdxmin2
subroutine cdxmin2(dtst, rhosig, iel, ifc, ifaxis)
Definition: subs1.f:751
get_strs_mask
subroutine get_strs_mask(mask, nxc, nzc, nel)
Definition: subs1.f:1757
fcaver
subroutine fcaver(xaver, a, iel, iface1)
Definition: subs1.f:870
axsf_e_2d
subroutine axsf_e_2d(au, av, u, v, h1, h2, ur, e)
Definition: subs1.f:2121
crs_strs
subroutine crs_strs(u1, u2, u3, v1, v2, v3)
Definition: subs1.f:1876
surf_mean
function surf_mean(u, ifld, bc_in, ierr)
Definition: subs1.f:1607
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342
unorm
subroutine unorm
Definition: subs1.f:352
axstrs
subroutine axstrs(a1, a2, a3, p1, p2, p3, h1, h2, matmod, nel)
Definition: subs1.f:1797
chktmg
subroutine chktmg(tol, res, w1, w2, mult, mask, imesh)
Definition: subs1.f:394
facdot
function facdot(A, B, IFACE1)
Definition: subs1.f:834
strs_orthok
subroutine strs_orthok(x, b, n, k, h1, h2, wt, ifld, w, ierr, matmod)
Definition: subs1.f:2513
axsf_fast
subroutine axsf_fast(au, av, aw, u, v, w, h1, h2, ifld)
Definition: subs1.f:2178
axstrs_nds
subroutine axstrs_nds(a1, a2, a3, p1, p2, p3, h1, h2, matmod, nel)
Definition: subs1.f:1807
hmhzsf
subroutine hmhzsf(name, u1, u2, u3, r1, r2, r3, h1, h2, rmask1, rmask2, rmask3, rmult, tol, maxit, matmod)
Definition: subs1.f:1096
strs_ortho_one
subroutine strs_ortho_one(x, b, n, k, h1, h2, wt, ifld, w, ierr, matmod)
Definition: subs1.f:2578
opnorm2w
function opnorm2w(v1, v2, v3, w)
Definition: subs1.f:2360
stnrate
subroutine stnrate(u1, u2, u3, nel, matmod)
Definition: subs1.f:1319
get_local_crs_galerkin_strs
subroutine get_local_crs_galerkin_strs(a, ncl, nxc, h1, h2, matmod)
Definition: subs1.f:1818
faccl3
subroutine faccl3(a, b, c, iface1)
Definition: subs1.f:944
aijuj
subroutine aijuj(au1, au2, au3, nel, ifaxis)
Definition: subs1.f:1435
uxyz
subroutine uxyz(u, ex, ey, ez, nel)
Definition: subs1.f:1458
strs_ortho_all
subroutine strs_ortho_all(x, b, n, k, h1, h2, wt, ifld, w, ierr, matmod)
Definition: subs1.f:2598
diagnos
subroutine diagnos
Definition: subs1.f:1078
set_vert_strs
subroutine set_vert_strs(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: subs1.f:1722
axiezz
subroutine axiezz(u2, eyy, ezz, nel)
Definition: subs1.f:1534
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
cdxmin3
subroutine cdxmin3(dtst, rhosig, iel, ifc)
Definition: subs1.f:797
setdtc
subroutine setdtc
Definition: subs1.f:463
ttxyz
subroutine ttxyz(ff, tx, ty, tz, nel)
Definition: subs1.f:2206
setsolv
subroutine setsolv
Definition: subs1.f:1083
stress
subroutine stress(h1, h2, nel, matmod, ifaxis)
Definition: subs1.f:1365
setprop
subroutine setprop
Definition: subs1.f:2618
fdm_h1a
subroutine fdm_h1a(z, r, d, nel, kt, rr)
Definition: subs1.f:1643
cggosf
subroutine cggosf(u1, u2, u3, r1, r2, r3, h1, h2, rmult, binv, vol, tin, maxit, matmod)
Definition: subs1.f:4
urst
subroutine urst(u, ur, us, ut, nel)
Definition: subs1.f:1490
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341
emerxit
subroutine emerxit
Definition: subs2.f:357
opadds
subroutine opadds(A1, A2, A3, B1, B2, B3, CONST, N, ISC)
Definition: subs2.f:101
op_glsc2_wt
function op_glsc2_wt(b1, b2, b3, x1, x2, x3, wt)
Definition: subs2.f:2258
facexv
subroutine facexv(A1, A2, A3, B1, B2, B3, IFACE1, IOP)
Definition: subs2.f:170
opdot
subroutine opdot(DP, A1, A2, A3, B1, B2, B3, N)
Definition: subs2.f:80
add3s
subroutine add3s(A, B, C, CONST, N)
Definition: subs2.f:349
op_vlsc2_wt
function op_vlsc2_wt(b1, b2, b3, x1, x2, x3, wt)
Definition: subs2.f:2234
setaxw1
subroutine setaxw1(IFAXWG)
Definition: subs2.f:2
facexs
subroutine facexs(A, B, IFACE1, IOP)
Definition: subs2.f:125
rmask
subroutine rmask(R1, R2, R3, NEL)
Definition: subs2.f:1801
vprops
subroutine vprops
Definition: vprops.f:2