KTH framework for Nek5000 toolboxes; testing version  0.0.1
induct.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2 c
3 c To do:
4 c
5 c Differing BC's imposed for ophinv, incomprn, etc.
6 c
7 c 1-shot Fast solver for Helmholtz and pressure
8 c
9 c
10 c-----------------------------------------------------------------------
11  subroutine induct (igeom)
12 c
13 c Solve the convection-diffusion equation for the B-field, with
14 c projection onto a div-free space.
15 c
16 c
17  include 'SIZE'
18  include 'INPUT'
19  include 'EIGEN'
20  include 'SOLN'
21  include 'TSTEP'
22  include 'MASS'
23 
24  common /scrns/ resv1(lx1,ly1,lz1,lelv)
25  $ , resv2(lx1,ly1,lz1,lelv)
26  $ , resv3(lx1,ly1,lz1,lelv)
27  $ , dv1(lx1,ly1,lz1,lelv)
28  $ , dv2(lx1,ly1,lz1,lelv)
29  $ , dv3(lx1,ly1,lz1,lelv)
30  common /scrvh/ h1(lx1,ly1,lz1,lelv)
31  $ , h2(lx1,ly1,lz1,lelv)
32 
33  ifield = ifldmhd
34 
35  if (igeom.eq.1) then ! old geometry, old velocity
36 
37  call makebsource_mhd
38 
39  else
40 
41  call lagbfield
42  call lagvel
43 
44  call elsasserh(igeom)
45 
46  call vol_flow ! check for fixed flow rate
47 
48 
49  endif
50 c
51  return
52  end
53 c--------------------------------------------------------------------
54  subroutine lagbfield
55 c
56 c Keep old B-field(s)
57 c
58  include 'SIZE'
59  include 'INPUT'
60  include 'SOLN'
61  include 'TSTEP'
62 c
63  do ilag=nbdinp-1,2,-1
64  call opcopy
65  $ ( bxlag(1,ilag ),bylag(1,ilag ),bzlag(1,ilag )
66  $ , bxlag(1,ilag-1),bylag(1,ilag-1),bzlag(1,ilag-1) )
67  enddo
68  call opcopy (bxlag,bylag,bzlag,bx,by,bz)
69 c
70  return
71  end
72 c--------------------------------------------------------------------
73  subroutine makebsource_mhd
74 c
75 c Make rhs for induction equation
76 c
77  include 'SIZE'
78  include 'SOLN'
79  include 'MASS'
80  include 'INPUT'
81  include 'TSTEP'
82  include 'CTIMER'
83 
84  if (icalld.eq.0) tbmhd=0.0
85  icalld = icalld+1
86  nbmhd = icalld
87  etime1 = dnekclock()
88 c
89  ifield = 1
90  call makeuf
91 c
92  ifield = ifldmhd
93  call makeufb
94  if (ifaxis) then
95 c do ifield = 2,3
96  do ifield = 2,npscal+1
97  call makeuq ! nonlinear terms
98  enddo
99  ifield = ifldmhd
100  endif
101 
102  if (ifnav.and.(.not.ifchar)) then
103  call advab_elsasser_fast
104  endif
105  if (ifchar) then
106  write(6,*) 'No IFCHAR for MHD, yet.'
107  call exitt
108  endif
109 
110  ifield = 1
111  if (iftran) call makeabf
112  call makebdf
113 c
114  ifield = ifldmhd
115  if (iftran) call makextb
116  call makebdfb
117 c
118  tbmhd=tbmhd+(dnekclock()-etime1)
119  return
120  end
121 c--------------------------------------------------------------------
122  subroutine makeufb
123 c
124 c Compute and add: (1) user specified forcing function (FX,FY,FZ)
125 c
126  include 'SIZE'
127  include 'SOLN'
128  include 'MASS'
129  include 'TSTEP'
130 C
131  time = time-dt
132  call nekuf (bmx,bmy,bmz)
133  call opcolv2 (bmx,bmy,bmz,vtrans(1,1,1,1,ifield),bm1)
134  time = time+dt
135 c
136  return
137  end
138 c--------------------------------------------------------------------
139  subroutine makextb
140 c
141 c Add extrapolation terms to magnetic source terms
142 c
143 c (nek5 equivalent for velocity is "makeabf")
144 c
145  include 'SIZE'
146  include 'INPUT'
147  include 'SOLN'
148  include 'MASS'
149  include 'TSTEP'
150 C
151  common /scrns/ ta1(lx1,ly1,lz1,lelv)
152  $ , ta2(lx1,ly1,lz1,lelv)
153  $ , ta3(lx1,ly1,lz1,lelv)
154 c
155  ntot1 = lx1*ly1*lz1*nelv
156 c
157  ab0 = ab(1)
158  ab1 = ab(2)
159  ab2 = ab(3)
160  call add3s2 (ta1,bbx1,bbx2,ab1,ab2,ntot1)
161  call add3s2 (ta2,bby1,bby2,ab1,ab2,ntot1)
162  call copy (bbx2,bbx1,ntot1)
163  call copy (bby2,bby1,ntot1)
164  call copy (bbx1,bmx,ntot1)
165  call copy (bby1,bmy,ntot1)
166  call add2s1 (bmx,ta1,ab0,ntot1)
167  call add2s1 (bmy,ta2,ab0,ntot1)
168  if (ldim.eq.3) then
169  call add3s2 (ta3,bbz1,bbz2,ab1,ab2,ntot1)
170  call copy (bbz2,bbz1,ntot1)
171  call copy (bbz1,bmz,ntot1)
172  call add2s1 (bmz,ta3,ab0,ntot1)
173  endif
174 c
175  return
176  end
177 c--------------------------------------------------------------------
178  subroutine makebdfb
179 C
180 C Add contributions to magnetic source from lagged BD terms.
181 C
182  include 'SIZE'
183  include 'SOLN'
184  include 'MASS'
185  include 'GEOM'
186  include 'INPUT'
187  include 'TSTEP'
188 C
189  COMMON /scrns/ ta1(lx1,ly1,lz1,lelv)
190  $ , ta2(lx1,ly1,lz1,lelv)
191  $ , ta3(lx1,ly1,lz1,lelv)
192  $ , tb1(lx1,ly1,lz1,lelv)
193  $ , tb2(lx1,ly1,lz1,lelv)
194  $ , tb3(lx1,ly1,lz1,lelv)
195  $ , h2(lx1,ly1,lz1,lelv)
196 C
197  ntot1 = lx1*ly1*lz1*nelv
198  const = 1./dt
199  CALL cmult2(h2,vtrans(1,1,1,1,ifield),const,ntot1)
200  CALL opcolv3c (tb1,tb2,tb3,bx,by,bz,bm1,bd(2))
201 C
202  DO 100 ilag=2,nbd
203  IF (ifgeom) THEN
204  CALL opcolv3c(ta1,ta2,ta3,bxlag(1,ilag-1),
205  $ bylag(1,ilag-1),
206  $ bzlag(1,ilag-1),
207  $ bm1lag(1,1,1,1,ilag-1),bd(ilag+1))
208  ELSE
209  CALL opcolv3c(ta1,ta2,ta3,bxlag(1,ilag-1),
210  $ bylag(1,ilag-1),
211  $ bzlag(1,ilag-1),
212  $ bm1 ,bd(ilag+1))
213  ENDIF
214  CALL opadd2 (tb1,tb2,tb3,ta1,ta2,ta3)
215  100 CONTINUE
216  CALL opadd2col (bmx,bmy,bmz,tb1,tb2,tb3,h2)
217 c
218  return
219  end
220 c--------------------------------------------------------------------
221  subroutine cresvib(resv1,resv2,resv3,h1,h2)
222 c
223 c Account for inhomogeneous Dirichlet boundary contributions
224 c in rhs of induction eqn.
225 c n
226 c Also, subtract off best estimate of grad p
227 c
228  include 'SIZE'
229  include 'TOTAL'
230  real resv1 (lx1,ly1,lz1,1)
231  real resv2 (lx1,ly1,lz1,1)
232  real resv3 (lx1,ly1,lz1,1)
233  real h1 (lx1,ly1,lz1,1)
234  real h2 (lx1,ly1,lz1,1)
235  common /scruz/ w1(lx1,ly1,lz1,lelv)
236  $ , w2(lx1,ly1,lz1,lelv)
237  $ , w3(lx1,ly1,lz1,lelv)
238 c
239 c
240  call bcdirvc (bx,by,bz,b1mask,b2mask,b3mask)
241  call extrapp (pm,pmlag)
242  call opgradt (resv1,resv2,resv3,pm)
243  call opadd2 (resv1,resv2,resv3,bmx,bmy,bmz)
244 c call opcopy (resv1,resv2,resv3,bmx,bmy,bmz)
245  call ophx (w1,w2,w3,bx,by,bz,h1,h2)
246  call opsub2 (resv1,resv2,resv3,w1,w2,w3)
247 c
248  return
249  end
250 c--------------------------------------------------------------------
251  subroutine cresvibp(resv1,resv2,resv3,h1,h2)
252 c
253 c Account for inhomogeneous Dirichlet boundary contributions
254 c in rhs of momentum eqn.
255 c n
256 c Also, subtract off best estimate of grad p
257 c
258  include 'SIZE'
259  include 'TOTAL'
260  real resv1 (lx1,ly1,lz1,1)
261  real resv2 (lx1,ly1,lz1,1)
262  real resv3 (lx1,ly1,lz1,1)
263  real h1 (lx1,ly1,lz1,1)
264  real h2 (lx1,ly1,lz1,1)
265  common /scruz/ w1(lx1,ly1,lz1,lelv)
266  $ , w2(lx1,ly1,lz1,lelv)
267  $ , w3(lx1,ly1,lz1,lelv)
268 c
269  call bcdirvc (vx,vy,vz,v1mask,v2mask,v3mask)
270  if (ifstrs) call bcneutr
271  call extrapp (pr,prlag)
272  call opgradt (resv1,resv2,resv3,pr)
273  call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
274 c call opcopy (resv1,resv2,resv3,bfx,bfy,bfz)
275  call ophx (w1,w2,w3,vx,vy,vz,h1,h2)
276  call opsub2 (resv1,resv2,resv3,w1,w2,w3)
277 c
278  return
279  end
280 c--------------------------------------------------------------------
281  subroutine incomprn (ux,uy,uz,up)
282 c
283 c Project U onto the closest incompressible field
284 c
285 c Input: U := (ux,uy,uz)
286 c
287 c Output: updated values of U, iproj, proj; and
288 c up := pressure currection req'd to impose div U = 0
289 c
290 c
291 c Dependencies: ifield ==> which "density" (vtrans) is used.
292 c
293 c Notes 1. up is _not_ scaled by bd(1)/dt. This should be done
294 c external to incompr().
295 c
296 c 2. up accounts _only_ for the perturbation pressure,
297 c not the current pressure derived from extrapolation.
298 c
299 c
300  include 'SIZE'
301  include 'TOTAL'
302  include 'CTIMER'
303 c
304  common /scrns/ w1(lx1,ly1,lz1,lelv)
305  $ , w2(lx1,ly1,lz1,lelv)
306  $ , w3(lx1,ly1,lz1,lelv)
307  $ , dv1(lx1,ly1,lz1,lelv)
308  $ , dv2(lx1,ly1,lz1,lelv)
309  $ , dv3(lx1,ly1,lz1,lelv)
310  $ , dp(lx2,ly2,lz2,lelv)
311  common /scrvh/ h1(lx1,ly1,lz1,lelv)
312  $ , h2(lx1,ly1,lz1,lelv)
313  common /scrhi/ h2inv(lx1,ly1,lz1,lelv)
314 
315  parameter(nset = 1 + lbelv/lelv)
316  common /orthov/ pset(lx2*ly2*lz2*lelv*mxprev,nset)
317  common /orthbi/ nprv(2)
318  logical ifprjp
319 
320  ifprjp=.false. ! Project out previous pressure solutions?
321  istart=param(95)
322  if (istep.ge.istart.and.istart.ne.0) ifprjp=.true.
323 
324  if (icalld.eq.0) tpres=0.0
325  icalld = icalld+1
326  npres = icalld
327  etime1 = dnekclock()
328 
329  ntot1 = lx1*ly1*lz1*nelv
330  ntot2 = lx2*ly2*lz2*nelv
331  intype = 1
332 
333  call rzero (h1,ntot1)
334  call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
335  call invers2 (h2inv,h2,ntot1)
336 
337  call opdiv (dp,ux,uy,uz)
338 
339  bdti = -bd(1)/dt
340  call cmult (dp,bdti,ntot2)
341 
342  call add2col2(dp,bm2,usrdiv,ntot2) ! User-defined divergence.
343 
344  call ortho (dp)
345 
346  i = 1 + ifield/ifldmhd
347  if (ifprjp) call setrhsp (dp,h1,h2,h2inv,pset(1,i),nprv(i))
348  scaledt = dt/bd(1)
349  scaledi = 1./scaledt
350  call cmult(dp,scaledt,ntot2) ! scale for tol
351  call esolver (dp,h1,h2,h2inv,intype)
352  call cmult(dp,scaledi,ntot2)
353  if (ifprjp) call gensolnp (dp,h1,h2,h2inv,pset(1,i),nprv(i))
354 
355  call add2(up,dp,ntot2)
356 
357  call opgradt (w1 ,w2 ,w3 ,dp)
358  call opbinv (dv1,dv2,dv3,w1 ,w2 ,w3 ,h2inv)
359  dtb = dt/bd(1)
360  call opadd2cm (ux ,uy ,uz ,dv1,dv2,dv3, dtb )
361 
362  if (ifmhd) call chkptol ! to avoid repetition
363 
364  tpres=tpres+(dnekclock()-etime1)
365 
366  return
367  end
368 c-----------------------------------------------------------------------
369  subroutine extrapp(p,plag)
370 C
371 C Pressure extrapolation
372 C
373  include 'SIZE'
374  include 'SOLN'
375  include 'TSTEP'
376 
377  real p (lx2,ly2,lz2,1)
378  $ ,plag (lx2,ly2,lz2,1)
379 
380  common /cgeom/ igeom
381 
382  ntot2 = lx2*ly2*lz2*nelv
383 
384  if (nbd.eq.2.and.nbdinp.gt.2.and.igeom.le.2) then
385  call copy(plag,p,ntot2)
386  elseif (nbd.eq.3.and.igeom.le.2) then
387 
388  const = dtlag(1)/dtlag(2)
389 
390  do i=1,ntot2
391  pnm1 = p(i,1,1,1)
392  pnm2 = plag(i,1,1,1)
393  p(i,1,1,1) = pnm1 + const*(pnm1-pnm2)
394  plag(i,1,1,1) = pnm1
395  enddo
396 
397  elseif (nbd.gt.3) then
398  WRITE (6,*) 'Pressure extrapolation cannot be completed'
399  WRITE (6,*) 'Try a lower-order temporal scheme'
400  call exitt
401  endif
402  return
403  end
404 c-----------------------------------------------------------------------
405  subroutine opzero(ux,uy,uz)
406  include 'SIZE'
407  include 'TOTAL'
408  real ux(1),uy(1),uz(1)
409 c
410  n = lx1*ly1*lz1*nelfld(ifield)
411  call rzero(ux,n)
412  call rzero(uy,n)
413  if (if3d) call rzero(uz,n)
414 c
415  return
416  end
417 c-----------------------------------------------------------------------
418  subroutine opnorm(unorm,ux,uy,uz,type3)
419  include 'SIZE'
420  include 'TOTAL'
421  character*3 type3
422  real ux(1),uy(1),uz(1)
423  real un(3),wn(3)
424 c
425  n = lx1*ly1*lz1*nelfld(ifield)
426  if (type3.eq.'L2 ') then
427  if (if3d) then
428  un(1) = vlsc3(ux,ux,bm1,n)
429  un(2) = vlsc3(uy,uy,bm1,n)
430  un(3) = vlsc3(uz,uz,bm1,n)
431  un(1) = un(1) + un(2) + un(3)
432  unorm = glsum(un(1),1)
433  if (unorm.gt.0) unorm = sqrt(unorm/volvm1)
434  else
435  un(1) = vlsc3(ux,ux,bm1,n)
436  un(2) = vlsc3(uy,uy,bm1,n)
437  un(1) = un(1) + un(2)
438  unorm = glsum(un(1),1)
439  if (unorm.gt.0) unorm = sqrt(unorm/volvm1)
440  endif
441  endif
442 c
443  return
444  end
445 c-----------------------------------------------------------------------
446  subroutine lorentz_force (lf,b1,b2,b3,w1,w2)
447 c
448 c Compute Lorentz force
449 c
450 c Input: B := (b1,b2,b3)
451 c
452 c Output: lf(1,ldim)
453 c
454 c Work arrays: w1(ltot) and w2(ltot)
455 c
456 c The output will not be continuous. However, it will be in
457 c the form appropriate for incorporation as a body force term
458 c in the variational formulation of the Navier-Stokes equations.
459 c
460 c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
461 c
462  include 'SIZE'
463 c
464  real lf(lx1*ly1*lz1*lelv,ldim)
465  real b1(lx1*ly1*lz1*lelv)
466  real b2(lx1*ly1*lz1*lelv)
467  real b3(lx1*ly1*lz1*lelv)
468 c
469  call curl(lf,b1,b2,b3,.false.,w1,w2)
470 c
471  ntot = lx1*ly1*lz1*nelv
472 c
473  do i=1,ntot
474  c1 = lf(i,2)*b3(i) - lf(i,3)*b2(i)
475  c2 = lf(i,3)*b1(i) - lf(i,1)*b3(i)
476  c3 = lf(i,1)*b2(i) - lf(i,2)*b1(i)
477  lf(i,1) = c1
478  lf(i,2) = c2
479  lf(i,3) = c3
480  enddo
481 c
482  return
483  end
484 c-----------------------------------------------------------------------
485  subroutine curl(vort,u,v,w,ifavg,work1,work2)
486 c
487  include 'SIZE'
488  include 'TOTAL'
489 c
490  logical ifavg
491 c
492  parameter(lt=lx1*ly1*lz1*lelv)
493  real vort(lt,3),work1(1),work2(1),u(1),v(1),w(1)
494 c
495  ntot = lx1*ly1*lz1*nelv
496  if (if3d) then
497 c work1=dw/dy ; work2=dv/dz
498  call dudxyz(work1,w,rym1,sym1,tym1,jacm1,1,2)
499  call dudxyz(work2,v,rzm1,szm1,tzm1,jacm1,1,3)
500  call sub3(vort(1,1),work1,work2,ntot)
501 c work1=du/dz ; work2=dw/dx
502  call dudxyz(work1,u,rzm1,szm1,tzm1,jacm1,1,3)
503  call dudxyz(work2,w,rxm1,sxm1,txm1,jacm1,1,1)
504  call sub3(vort(1,2),work1,work2,ntot)
505 c work1=dv/dx ; work2=du/dy
506  call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
507  call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
508  call sub3(vort(1,3),work1,work2,ntot)
509  else
510 c work1=dv/dx ; work2=du/dy
511  call dudxyz(work1,v,rxm1,sxm1,txm1,jacm1,1,1)
512  call dudxyz(work2,u,rym1,sym1,tym1,jacm1,1,2)
513  call sub3(vort(1,3),work1,work2,ntot)
514  endif
515 c
516 c Avg at bndry
517 c
518  if (ifavg) then
519  ifielt = ifield
520  ifield = 1
521  if (if3d) then
522  do idim=1,ldim
523  call col2 (vort(1,idim),bm1,ntot)
524  call dssum (vort(1,idim),lx1,ly1,lz1)
525  call col2 (vort(1,idim),binvm1,ntot)
526  enddo
527  else
528  call col2 (vort(1,3),bm1,ntot) ! NOTE: This differs from
529  call dssum (vort(1,3),lx1,ly1,lz1) ! "comp_vort", which returns
530  call col2 (vort(1,3),binvm1,ntot) ! vorticity as 1st entry in vort
531  endif
532  ifield = ifielt
533  endif
534 c
535  return
536  end
537 c-----------------------------------------------------------------------
538  subroutine lorentz_force2(lf,b1,b2,b3)
539 c
540 c Compute Lorentz force
541 c
542 c Input: B := (b1,b2,b3)
543 c
544 c Output: lf(1,ldim)
545 c
546 c Work arrays: w1(ltot) and w2(ltot)
547 c
548 c The output will not be continuous. However, it will be in
549 c the form appropriate for incorporation as a body force term
550 c in the variational formulation of the Navier-Stokes equations.
551 c
552 c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
553 c
554  include 'SIZE'
555 c
556  real lf(lx1*ly1*lz1*ldim,lelt)
557  real b1(lx1*ly1*lz1,lelt)
558  real b2(lx1*ly1*lz1,lelt)
559  real b3(lx1*ly1*lz1,lelt)
560 c
561  integer e
562 c
563  do e = 1,nelt ! NOTE: the order is different from v. 1
564  call lorentz_force_e(lf(1,e),b1(1,e),b2(1,e),b3(1,e),e)
565  enddo
566 c
567  return
568  end
569 c-----------------------------------------------------------------------
570  subroutine lorentz_force_e(lf,b1,b2,b3,e)
571 c
572 c Compute Lorentz force field for a single element
573 c
574 c Input: B := (b1,b2,b3)
575 c
576 c Output: lf(1,ldim)
577 c
578 c Work arrays: cb(lxyzd) and w2(lxyzd)
579 c
580 c The output will not be continuous. However, it will be in
581 c the form appropriate for incorporation as a body force term
582 c in the variational formulation of the Navier-Stokes equations.
583 c (i.e., rhs(NS) = rhs(NS) + B*lf, where B is the mass matrix)
584 c
585 c Dealiasing strategy:
586 c
587 c e e -1 ~ T ~e ~ ~
588 c lf = ( B ) I B ( I j x I b )
589 c i ~ i
590 c
591 c ~
592 c Here, I is the interpolant from N to M, where M = 3/2 N.
593 c
594 c ~ -1
595 c j is a special curl (sans Jacobian )
596 c
597 c b is the B-field on the N-points
598 c
599 c ~e
600 c B is the local mass matrix on the M-points (sans Jacabian)
601 c
602 c e
603 c B is the local mass matrix on the N-points (with Jacabian)
604 c
605 c The last B is req'd to compensate for the subsequent multiplication
606 c by B that takes place once the rhs is formed.
607 c
608 c
609 c
610  include 'SIZE'
611  include 'DEALIAS'
612  include 'GEOM'
613  include 'WZ'
614 c
615  real lf(lx1*ly1*lz1,3)
616  real b1(lx1*ly1*lz1)
617  real b2(lx1*ly1*lz1)
618  real b3(lx1*ly1*lz1)
619  integer d,e
620 c
621  integer icalld
622  save icalld
623  data icalld /0/
624 c
625  common /ctmp1x/ lfd(lxd*lyd*lzd,3)
626  $ , bd(lxd*lyd*lzd,3)
627  $ , cb(lx1*ly1*lz1,3)
628  $ , cbd(lxd*lyd*lzd,3)
629  real lfd
630 
631  if (icalld .eq. 0) then
632  write(6,*) 'CALL SET PROJ',lx1,lxd
633  call setmap (lx1,lxd) ! Set up interpolation operators
634  call setproj(lx1,lxd) ! Set up interpolation operators
635  icalld = icalld + 1
636  endif
637 c
638  call spec_curl_e (cb,b1,b2,b3 !local curl, w/o Jacobian
639  $ , rxm1(1,1,1,e),rym1(1,1,1,e),rzm1(1,1,1,e)
640  $ , sxm1(1,1,1,e),sym1(1,1,1,e),szm1(1,1,1,e)
641  $ , txm1(1,1,1,e),tym1(1,1,1,e),tzm1(1,1,1,e) )
642 c
643  do d=1,ldim ! interpolate to M points
644  call specmp(cbd(1,d),lxd,cb(1,d),lx1,im1d,im1dt,bd)
645  enddo
646  call specmp(bd(1,1),lxd,b1,lx1,im1d,im1dt,lfd)
647  call specmp(bd(1,2),lxd,b2,lx1,im1d,im1dt,lfd)
648  call specmp(bd(1,3),lxd,b3,lx1,im1d,im1dt,lfd)
649 c
650  nxyzd = lxd*lyd*lzd
651  do i=1,nxyzd
652  lfd(i,1) = cbd(i,2)*bd(i,3) - cbd(i,3)*bd(i,2) ! Curl B x B
653  lfd(i,2) = cbd(i,3)*bd(i,1) - cbd(i,1)*bd(i,3)
654  lfd(i,3) = cbd(i,1)*bd(i,2) - cbd(i,2)*bd(i,1)
655  enddo
656 c
657 c Project back and simultaneous collocate with local quadrature weights
658 c
659 c ~ ~ ~
660 c P := Pz x Py x Px = Iz*B x Iy*B x Ix*B
661 c
662 c
663 c ~ M
664 c where B = diag (rho )
665 c i
666 c
667  call specmp(lf(1,1),lx1,lfd(1,1),lxd,pmd1,pmd1t,cbd)
668  call specmp(lf(1,2),lx1,lfd(1,2),lxd,pmd1,pmd1t,cbd)
669  call specmp(lf(1,3),lx1,lfd(1,3),lxd,pmd1,pmd1t,cbd)
670 c
671 c Finally, divide by local mass matrix in anticipation of subsequent
672 c multiply by BM1.
673 c
674  nxyz = lx1*ly1*lz1
675  do i=1,nxyz
676 c scale = 1./(w3m1(i,1,1)*jacm1(i,1,1,e))
677  scale = 1./jacm1(i,1,1,e)
678  lf(i,1) = scale*lf(i,1)
679  lf(i,2) = scale*lf(i,2)
680  lf(i,3) = scale*lf(i,3)
681  enddo
682 c
683  return
684  end
685 c-----------------------------------------------------------------------
686  subroutine spec_curl_e (cb,b1,b2,b3,rx,ry,rz,sx,sy,sz,tx,ty,tz)
687 c
688 c local curl, multiplied by Jacobian
689 c
690  include 'SIZE'
691  include 'DXYZ'
692 c
693  real cb(lx1*ly1*lz1,3) ! Output J*curl B (J:=Jacobian)
694  real b1(1),b2(1),b3(1) ! Input B-field
695 c
696  real rx(1),ry(1),rz(1) ! Metrics
697  real sx(1),sy(1),sz(1)
698  real tx(1),ty(1),tz(1)
699 c
700  common /ctmp0x/ br(lx1*ly1*lz1),bs(lx1*ly1*lz1),bt(lx1*ly1*lz1)
701 c
702 c / db3 db2 \
703 c cb1 = J | --- - --- | ! Keep J ( J:=Jacobian)
704 c \ dy dz /
705 c
706 c
707 c / db1 db3 \
708 c cb2 = J | --- - --- | ! Keep J ( J:=Jacobian)
709 c \ dz dx /
710 c
711 c
712 c / db2 db1 \
713 c cb3 = J | --- - --- | ! Keep J ( J:=Jacobian)
714 c \ dx dy /
715 c
716 c
717 c Note:
718 c
719 c db2 db2 dr db2 ds db2 dt
720 c J --- = --- J -- + --- J -- + --- J --
721 c dz dr dz ds dz dt dz
722 c
723 c etc.
724 c
725 c
726 c
727  nxyz = lx1*ly1*lz1
728 c
729  n=lx1-1
730  call local_grad3(br,bs,bt,b1,n,1,dxm1,dxtm1)
731  do i=1,nxyz
732  cb(i,2) = (br(i)*rz(i)+bs(i)*sz(i)+bt(i)*tz(i))
733  cb(i,3) = -(br(i)*ry(i)+bs(i)*sy(i)+bt(i)*ty(i))
734  enddo
735 c
736  call local_grad3(br,bs,bt,b2,n,1,dxm1,dxtm1)
737  do i=1,nxyz
738  cb(i,1) = -(br(i)*rz(i)+bs(i)*sz(i)+bt(i)*tz(i))
739  cb(i,3) = (br(i)*rx(i)+bs(i)*sx(i)+bt(i)*tx(i))
740  $ + cb(i,3)
741  enddo
742 c
743  call local_grad3(br,bs,bt,b3,n,1,dxm1,dxtm1)
744  do i=1,nxyz
745  cb(i,1) = (br(i)*ry(i)+bs(i)*sy(i)+bt(i)*ty(i))
746  $ + cb(i,1)
747  cb(i,2) = -(br(i)*rx(i)+bs(i)*sx(i)+bt(i)*tx(i))
748  $ + cb(i,2)
749  enddo
750 c
751  return
752  end
753 c-----------------------------------------------------------------------
754  subroutine specx(b,nb,a,na,ba,ab,w)
755 c
756  include 'SIZE'
757  include 'INPUT'
758  real b(1),a(1)
759  real w(1)
760 c
761  n=na*na*na
762  do i=1,n
763  b(i) = a(i)
764  enddo
765 c
766  return
767  end
768 c-----------------------------------------------------------------------
769  subroutine phys_to_elsasser(u1,u2,u3,b1,b2,b3,n)
770 c
771  real u1(1),u2(1),u3(1),b1(1),b2(1),b3(1)
772 c
773  do i=1,n
774 c
775  zpx = u1(i) + b1(i)
776  zpy = u2(i) + b2(i)
777  zpz = u3(i) + b3(i)
778 c
779  zmx = u1(i) - b1(i)
780  zmy = u2(i) - b2(i)
781  zmz = u3(i) - b3(i)
782 c
783  u1(i) = zpx
784  u2(i) = zpy
785  u3(i) = zpz
786 c
787  b1(i) = zmx
788  b2(i) = zmy
789  b3(i) = zmz
790 c
791  enddo
792 c
793  return
794  end
795 c-----------------------------------------------------------------------
796  subroutine elsasser_to_phys(u1,u2,u3,b1,b2,b3,n)
797 c
798  real u1(1),u2(1),u3(1),b1(1),b2(1),b3(1)
799 c
800  do i=1,n
801 c
802  zpx = 0.5*( u1(i) + b1(i) )
803  zpy = 0.5*( u2(i) + b2(i) )
804  zpz = 0.5*( u3(i) + b3(i) )
805 c
806  zmx = 0.5*( u1(i) - b1(i) )
807  zmy = 0.5*( u2(i) - b2(i) )
808  zmz = 0.5*( u3(i) - b3(i) )
809 c
810  u1(i) = zpx
811  u2(i) = zpy
812  u3(i) = zpz
813 c
814  b1(i) = zmx
815  b2(i) = zmy
816  b3(i) = zmz
817 c
818  enddo
819 c
820  return
821  end
822 c-----------------------------------------------------------------------
823  subroutine phys_to_elsasser2(u1,b1,n)
824 c
825  real u1(1),b1(1)
826 c
827  do i=1,n
828  zpx = u1(i) + b1(i)
829  zmx = u1(i) - b1(i)
830  u1(i) = zpx
831  b1(i) = zmx
832  enddo
833 c
834  return
835  end
836 c-----------------------------------------------------------------------
837  subroutine elsasser_to_phys2(u1,b1,n)
838 c
839  real u1(1),b1(1)
840 c
841  do i=1,n
842  zpx = 0.5*( u1(i) + b1(i) )
843  zmx = 0.5*( u1(i) - b1(i) )
844  u1(i) = zpx
845  b1(i) = zmx
846  enddo
847 c
848  return
849  end
850 c-----------------------------------------------------------------------
851  subroutine elsasserh(igeom)
852 c
853 c
854 c Solve MHD in Elsasser variables
855 c
856 c
857  include 'SIZE'
858  include 'INPUT'
859  include 'EIGEN'
860  include 'SOLN'
861  include 'TSTEP'
862  include 'MASS'
863  include 'GEOM'
864 C
865  common /scrnt/ besv1(lbx1,lby1,lbz1,lbelv)
866  $ , besv2(lbx1,lby1,lbz1,lbelv)
867  $ , besv3(lbx1,lby1,lbz1,lbelv)
868  COMMON /scrns/ resv1(lx1,ly1,lz1,lelv)
869  $ , resv2(lx1,ly1,lz1,lelv)
870  $ , resv3(lx1,ly1,lz1,lelv)
871  $ , dv1(lx1,ly1,lz1,lelv)
872  $ , dv2(lx1,ly1,lz1,lelv)
873  $ , dv3(lx1,ly1,lz1,lelv)
874  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
875  $ , h2(lx1,ly1,lz1,lelv)
876 c
877  n = lx1*ly1*lz1*nelv
878 c
879 c New geometry, new velocity
880 c
881  intype = -1
882 c
883  ifield = 1
884  call sethlm (h1,h2,intype)
885  call cresvibp (resv1,resv2,resv3,h1,h2)
886 c
887  ifield = ifldmhd
888  call sethlm (h1,h2,intype)
889  call cresvib (besv1,besv2,besv3,h1,h2)
890 c
891 c
892  ifield = 1
893  call sethlm (h1,h2,intype)
894 
895  call ophinv (dv1,dv2,dv3,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
896 
897  call opadd2 (vx,vy,vz,dv1,dv2,dv3)
898 
899  if (filtertype.eq.1) then
900  alpha_filt=param(103) ! Optional Filtering
901  call q_filter(alpha_filt)
902  endif
903 
904  call incomprn (vx,vy,vz,pr) ! project U onto div-free space
905 
906  ifield = ifldmhd
907  call sethlm (h1,h2,intype)
908 
909  call ophinv (dv1,dv2,dv3,besv1,besv2,besv3,h1,h2,tolhv,nmxv)
910  call opadd2 (bx,by,bz,dv1,dv2,dv3)
911 
912  call incomprn (bx,by,bz,pm) ! project B onto div-free space
913 
914  return
915  end
916 c--------------------------------------------------------------------
917  subroutine compute_cfl(cfl,u,v,w,dt)
918 c
919 c Given velocity field (u,v,w) and dt, compute current CFL number.
920 c
921  include 'SIZE'
922  include 'GEOM'
923  include 'INPUT'
924  include 'WZ'
925  include 'SOLN'
926 c
927  real u(lx1,ly1,lz1,nelv),v(lx1,ly1,lz1,nelv),w(lx1,ly1,lz1,nelv)
928 c
929 c Store the inverse jacobian to speed up this operation
930 c
931  common /cfldx/ dri(lx1),dsi(ly1),dti(lz1)
932 c
933  integer e
934 c
935  integer icalld
936  save icalld
937  data icalld /0/
938 c
939  if (icalld.eq.0) then
940  icalld=1
941  call getdr(dri,zgm1(1,1),lx1)
942  call getdr(dsi,zgm1(1,2),ly1)
943  if (if3d) call getdr(dti,zgm1(1,3),lz1)
944  endif
945 
946  cfl = 0.
947  l = 0
948 
949  if (if3d) then
950  nxyz = lx1*ly1*lz1
951  do e=1,nelv
952  do k=1,lz1
953  do j=1,ly1
954  do i=1,lx1
955  l = l+1
956  ur = ( u(i,j,k,e)*rxm1(i,j,k,e)
957  $ + v(i,j,k,e)*rym1(i,j,k,e)
958  $ + w(i,j,k,e)*rzm1(i,j,k,e) ) * jacmi(l,1)
959  us = ( u(i,j,k,e)*sxm1(i,j,k,e)
960  $ + v(i,j,k,e)*sym1(i,j,k,e)
961  $ + w(i,j,k,e)*szm1(i,j,k,e) ) * jacmi(l,1)
962  ut = ( u(i,j,k,e)*txm1(i,j,k,e)
963  $ + v(i,j,k,e)*tym1(i,j,k,e)
964  $ + w(i,j,k,e)*tzm1(i,j,k,e) ) * jacmi(l,1)
965 
966  cflr = abs(dt*ur*dri(i))
967  cfls = abs(dt*us*dsi(j))
968  cflt = abs(dt*ut*dti(k))
969 
970  cflm = cflr + cfls + cflt
971  cfl = max(cfl,cflm)
972 
973  cflf(i,j,k,e) = cflm
974 
975  enddo
976  enddo
977  enddo
978  enddo
979  else
980  nxyz = lx1*ly1
981  do e=1,nelv
982  do j=1,ly1
983  do i=1,lx1
984  l = l+1
985  ur = ( u(i,j,1,e)*rxm1(i,j,1,e)
986  $ + v(i,j,1,e)*rym1(i,j,1,e) ) * jacmi(l,1)
987  us = ( u(i,j,1,e)*sxm1(i,j,1,e)
988  $ + v(i,j,1,e)*sym1(i,j,1,e) ) * jacmi(l,1)
989 
990  cflr = abs(dt*ur*dri(i))
991  cfls = abs(dt*us*dsi(j))
992 
993  cflm = cflr + cfls
994  cfl = max(cfl,cflm)
995 
996  cflf(i,j,1,e) = cflm
997 
998  enddo
999  enddo
1000  enddo
1001  endif
1002 c
1003  cfl = glmax(cfl,1)
1004 c
1005  return
1006  end
1007 c-----------------------------------------------------------------------
1008  subroutine getdr(dri,zgm1,lx1)
1009  real dri(lx1),zgm1(lx1)
1010 c
1011  dri(1) = zgm1(2) - zgm1(1) ! Compute 1/Dx
1012  do i=2,lx1-1
1013  dri(i) = 0.5*( zgm1(i+1) - zgm1(i-1) )
1014  enddo
1015  dri(lx1) = zgm1(lx1) - zgm1(lx1-1)
1016 c
1017  call invcol1(dri,lx1)
1018 c
1019  return
1020  end
1021 c-----------------------------------------------------------------------
1022  subroutine ophinv(o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
1023 C
1024 C Ok = (H1*A+H2*B)-1 * Ik (implicit)
1025 C
1026  include 'SIZE'
1027  include 'INPUT'
1028  include 'MASS'
1029  include 'SOLN'
1030  include 'VPROJ'
1031  include 'TSTEP'
1032 c logical ifproj
1033 
1034  real o1 (lx1,ly1,lz1,1) , o2 (lx1,ly1,lz1,1) , o3 (lx1,ly1,lz1,1)
1035  real i1 (lx1,ly1,lz1,1) , i2 (lx1,ly1,lz1,1) , i3 (lx1,ly1,lz1,1)
1036  real h1 (lx1,ly1,lz1,1) , h2 (lx1,ly1,lz1,1)
1037 
1038 c ifproj = .false.
1039 c if (param(94).gt.0) ifproj = .true.
1040 c if (ifprojfld(ifield)) ifproj = .true.
1041 c
1042 c if (.not.ifproj) then
1043 c if (ifield.eq.1) call ophinv
1044 c $ (o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
1045 c if (ifield.eq.ifldmhd) call ophinv
1046 c $ (o1,o2,o3,i1,i2,i3,h1,h2,tolh,nmxhi)
1047 c return
1048 c endif
1049 
1050  mtmp = param(93)
1051  do i=1,2*ldim
1052  ivproj(1,i) = min(mxprev,mtmp) - 1
1053  enddo
1054 
1055  imesh = 1
1056 
1057  if (ifstrs) then
1058  matmod = 0
1059  call hmhzsf ('NOMG',o1,o2,o3,i1,i2,i3,h1,h2,
1060  $ v1mask,v2mask,v3mask,vmult,
1061  $ tolh,nmxhi,matmod)
1062  else
1063  if (ifield.eq.1) then
1064  call hsolve ('VELX',o1,i1,h1,h2,v1mask,vmult
1065  $ ,imesh,tolh,nmxhi,1
1066  $ ,vproj(1,1),ivproj(1,1),binvm1)
1067  call hsolve ('VELY',o2,i2,h1,h2,v2mask,vmult
1068  $ ,imesh,tolh,nmxhi,2
1069  $ ,vproj(1,2),ivproj(1,2),binvm1)
1070  if (if3d)
1071  $ call hsolve ('VELZ',o3,i3,h1,h2,v3mask,vmult
1072  $ ,imesh,tolh,nmxhi,3
1073  $ ,vproj(1,3),ivproj(1,3),binvm1)
1074  elseif (ifield.eq.ifldmhd) then ! B-field
1075  call hsolve (' BX ',o1,i1,h1,h2,b1mask,vmult
1076  $ ,imesh,tolh,nmxhi,1
1077  $ ,vproj(1,4),ivproj(1,4),binvm1)
1078  call hsolve (' BY ',o2,i2,h1,h2,b2mask,vmult
1079  $ ,imesh,tolh,nmxhi,2
1080  $ ,vproj(1,5),ivproj(1,5),binvm1)
1081  if (if3d)
1082  $ call hsolve (' BZ ',o3,i3,h1,h2,b3mask,vmult
1083  $ ,imesh,tolh,nmxhi,3
1084  $ ,vproj(1,6),ivproj(1,6),binvm1)
1085  endif
1086  endif
1087 C
1088  return
1089  end
1090 c--------------------------------------------------------------------
1091  subroutine set_ifbcor
1092  include 'SIZE'
1093  include 'GEOM'
1094  include 'INPUT'
1095 c include 'TSTEP' ! ifield?
1096 
1097  common /nekcb/ cb
1098  character cb*3
1099 
1100  ifbcor = .true.
1101  ifield = ifldmhd
1102 
1103  nface = 2*ldim
1104  do iel=1,nelv
1105  do ifc=1,nface
1106  cb = cbc(ifc,iel,ifield)
1107  if (cb.eq.'ndd' .or. cb.eq.'dnd' .or. cb.eq.'ddn')
1108  $ ifbcor = .false.
1109  enddo
1110  enddo
1111 
1112  call gllog(ifbcor , .false.)
1113 
1114  if (nio.eq.0) write (6,*) 'IFBCOR =',ifbcor
1115 
1116  return
1117  end
1118 c--------------------------------------------------------------------
1119 c subroutine set_ifbcor_old(ifbcor)
1120 cc
1121 cc This is a quick hack for the rings problem - it is not general,
1122 cc but will also work fine for the periodic in Z problem
1123 cc
1124 c include 'SIZE'
1125 c include 'TOTAL'
1126 c
1127 c logical ifbcor
1128 c
1129 c integer e,f
1130 c character*1 cb1(3)
1131 c
1132 c itest = 0
1133 c
1134 c do e=1,nelv
1135 c do f=1,2*ldim
1136 c call chcopy(cb1,cbc(f,e,ifldmhd),3)
1137 c if (cb1(3).eq.'n'.or.cb1(3).eq.'N') itest = 1
1138 c enddo
1139 c enddo
1140 c
1141 c itest = iglmax(itest,1)
1142 c
1143 c ifbcor = .true. ! adjust mean pressure to rmv hydrostatic mode
1144 c
1145 c
1146 c if (itest.eq.1) ifbcor = .false. ! do not adjust mean pressure
1147 c
1148 c return
1149 c end
1150 c--------------------------------------------------------------------
1151  subroutine setrhsp(p,h1,h2,h2inv,pset,niprev)
1152 C
1153 C Project soln onto best fit in the "E" norm.
1154 C
1155  include 'SIZE'
1156  include 'INPUT'
1157  include 'MASS'
1158  include 'SOLN'
1159  include 'TSTEP'
1160 
1161  real p (lx2,ly2,lz2,lelv)
1162  real h1 (lx1,ly1,lz1,lelv)
1163  real h2 (lx1,ly1,lz1,lelv)
1164  real h2inv(lx1,ly1,lz1,lelv)
1165  real pset (lx2*ly2*lz2*lelv,mxprev)
1166 
1167  parameter(ltot2=lx2*ly2*lz2*lelv)
1168  common /orthox/ pbar(ltot2),pnew(ltot2)
1169  common /orthos/ alpha(mxprev),work(mxprev)
1170  common /orthoi/ nprev,mprev
1171 
1172  nprev = niprev
1173  if (nprev.eq.0) return
1174 
1175 c Diag to see how much reduction in the residual is attained.
1176  ntot2 = lx2*ly2*lz2*nelv
1177  alpha1 = glsc3(p,p,bm2inv,ntot2)
1178  if (alpha1.gt.0) then
1179  alpha1 = sqrt(alpha1/volvm2)
1180  else
1181  return
1182  endif
1183 
1184  call updrhse(p,h1,h2,h2inv,ierr) ! Update rhs's if E-matrix has changed
1185 c if (ierr.eq.1) Nprev=0 ! Doesn't happen w/ new formulation
1186 
1187  do i=1,nprev ! Perform Gram-Schmidt for previous soln's.
1188  alpha(i) = vlsc2(p,pset(1,i),ntot2)
1189  enddo
1190  call gop(alpha,work,'+ ',nprev)
1191 
1192  call rzero(pbar,ntot2)
1193  do i=1,nprev
1194  call add2s2(pbar,pset(1,i),alpha(i),ntot2)
1195  enddo
1196 C
1197  intetype = 1
1198  call cdabdtp(pnew,pbar,h1,h2,h2inv,intetype)
1199  call sub2 (p,pnew,ntot2)
1200 
1201 c ................................................................
1202  alpha2 = glsc3(p,p,bm2inv,ntot2) ! Diagnostics
1203  if (alpha2.gt.0) then
1204  alpha2 = sqrt(alpha2/volvm2)
1205  ratio = alpha1/alpha2
1206  n10=min(10,nprev)
1207 c if (nio.eq.0) write(6,11) istep,nprev,(alpha(i),i=1,n10)
1208 c if (nio.eq.0) write(6,12) istep,nprev,alpha1,alpha2,ratio
1209  11 format(2i5,' alpha:',1p10e12.4)
1210  12 format(i6,i4,1p3e12.4,' alph12')
1211 
1212  if (nio.eq.0) write(6,13) istep,' Project PRES ',
1213  & alpha2,alpha1,ratio,nprev,mxprev
1214  13 format(i11,a,6x,1p3e13.4,i4,i4)
1215  endif
1216 c ................................................................
1217 
1218  return
1219  end
1220 c-----------------------------------------------------------------------
1221  subroutine gensolnp(p,h1,h2,h2inv,pset,nprev)
1222 C
1223 C Reconstruct the solution to the original problem by adding back
1224 C the previous solutions
1225 C
1226  include 'SIZE'
1227  include 'INPUT'
1228 
1229  real p (lx2,ly2,lz2,lelv)
1230  real h1 (lx1,ly1,lz1,lelv)
1231  real h2 (lx1,ly1,lz1,lelv)
1232  real h2inv(lx1,ly1,lz1,lelv)
1233  real pset (lx2*ly2*lz2*lelv,mxprev)
1234 
1235  parameter(ltot2=lx2*ly2*lz2*lelv)
1236  common /orthox/ pbar(ltot2),pnew(ltot2)
1237  common /orthos/ alpha(mxprev),work(mxprev)
1238 
1239  mprev=param(93)
1240  mprev=min(mprev,mxprev)
1241 
1242  ntot2=lx2*ly2*lz2*nelv
1243 
1244  if (nprev.lt.mprev) then
1245  nprev = nprev+1
1246  call copy (pset(1,nprev),p,ntot2) ! Save current solution
1247  call add2 (p,pbar,ntot2) ! Reconstruct solution.
1248  call econjp(pset,nprev,h1,h2,h2inv,ierr) ! Orthonormalize set
1249 
1250  if (ierr.eq.1) then
1251  nprev = 1
1252  call copy (pset(1,nprev),p,ntot2) ! Save current solution
1253  call econjp(pset,nprev,h1,h2,h2inv,ierr) ! and orthonormalize.
1254  endif
1255 
1256  else ! (uses pnew).
1257  nprev = 1
1258  call add2 (p,pbar,ntot2) ! Reconstruct solution.
1259  call copy (pset(1,nprev),p,ntot2) ! Save current solution
1260  call econjp(pset,nprev,h1,h2,h2inv,ierr) ! and orthonormalize.
1261  endif
1262 
1263  return
1264  end
1265 c-----------------------------------------------------------------------
1266  subroutine econjp(pset,nprev,h1,h2,h2inv,ierr)
1267 
1268 c Orthogonalize the soln wrt previous soln's for which we already
1269 c know the soln.
1270 
1271  include 'SIZE'
1272  include 'INPUT'
1273 
1274  real p (lx2,ly2,lz2,lelv)
1275  real h1 (lx1,ly1,lz1,lelv)
1276  real h2 (lx1,ly1,lz1,lelv)
1277  real h2inv(lx1,ly1,lz1,lelv)
1278  real pset (lx2*ly2*lz2*lelv,mxprev)
1279 
1280  parameter(ltot2=lx2*ly2*lz2*lelv)
1281  common /orthox/ pbar(ltot2),pnew(ltot2)
1282  common /orthos/ alpha(mxprev),work(mxprev)
1283 
1284  ierr = 0
1285 
1286  ntot2=lx2*ly2*lz2*nelv
1287 
1288 C
1289 C Gram Schmidt, w re-orthogonalization
1290 C
1291  npass=1
1292 c if (abs(param(105)).eq.2) npass=2
1293  do ipass=1,npass
1294 
1295  intetype=1
1296  call cdabdtp(pnew,pset(1,nprev),h1,h2,h2inv,intetype)
1297  alphad = glsc2(pnew,pset(1,nprev),ntot2) ! compute part of the norm
1298 
1299  nprev1 = nprev-1
1300  do i=1,nprev1 ! Gram-Schmidt
1301  alpha(i) = vlsc2(pnew,pset(1,i),ntot2)
1302  enddo
1303  if (nprev1.gt.0) call gop(alpha,work,'+ ',nprev1)
1304 
1305  do i=1,nprev1
1306  alpham = -alpha(i)
1307  call add2s2(pset(1,nprev),pset(1,i),alpham,ntot2)
1308  alphad = alphad - alpha(i)**2
1309  enddo
1310  enddo
1311 C
1312 C .Normalize new element in P~
1313 C
1314  if (alphad.le.0) then
1315  write(6,*) .le.'ERROR: alphad 0 in econjp',alphad,nprev
1316  ierr = 1
1317  return
1318  endif
1319  alphad = 1./sqrt(alphad)
1320  call cmult(pset(1,nprev),alphad,ntot2)
1321 
1322  return
1323  end
1324 c-----------------------------------------------------------------------
1325  subroutine advab_elsasser_fast
1326 C
1327 C Eulerian scheme, add convection term to forcing function
1328 C at current time step.
1329 C
1330  include 'SIZE'
1331  include 'INPUT'
1332  include 'GEOM'
1333  include 'SOLN'
1334  include 'MASS'
1335  include 'TSTEP'
1336 c
1337  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
1338  common /scrns/ wk(2*ltd)
1339  $ , fx(lxy),fy(lxy),fz(lxy)
1340  $ , gx(lxy),gy(lxy),gz(lxy)
1341  $ , zr(ltd),zs(ltd),zt(ltd)
1342  $ , tr(ltd,3),zp(ltd,3),zm(ltd,3)
1343 
1344  integer e
1345  integer icalld
1346  save icalld
1347  data icalld /0/
1348 
1349  if (icalld.eq.0) call set_dealias_rx
1350  icalld = icalld + 1
1351 
1352  nxyz1 = lx1*ly1*lz1
1353  nxyzd = lxd*lyd*lzd
1354 
1355 
1356  if (if3d) then
1357 c
1358  do e=1,nelv
1359 
1360 c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
1361 
1362 c write(6,*) nid,' inside fast',e,nxyz1
1363 
1364  call add3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
1365  call intp_rstd(zp(1,1),wk,lx1,lxd,if3d,0) ! 0 --> forward
1366 
1367  call add3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
1368  call intp_rstd(zp(1,2),wk,lx1,lxd,if3d,0)
1369 
1370  call add3(wk,vz(1,1,1,e),bz(1,1,1,e),nxyz1)
1371  call intp_rstd(zp(1,3),wk,lx1,lxd,if3d,0)
1372 
1373  call sub3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
1374  call intp_rstd(zm(1,1),wk,lx1,lxd,if3d,0)
1375 
1376  call sub3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
1377  call intp_rstd(zm(1,2),wk,lx1,lxd,if3d,0)
1378 
1379  call sub3(wk,vz(1,1,1,e),bz(1,1,1,e),nxyz1)
1380  call intp_rstd(zm(1,3),wk,lx1,lxd,if3d,0)
1381 
1382  do i=1,nxyzd ! Convert convector (zm) to r-s-t coordinates
1383  tr(i,1)=
1384  $ rx(i,1,e)*zm(i,1)+rx(i,2,e)*zm(i,2)+rx(i,3,e)*zm(i,3)
1385  tr(i,2)=
1386  $ rx(i,4,e)*zm(i,1)+rx(i,5,e)*zm(i,2)+rx(i,6,e)*zm(i,3)
1387  tr(i,3)=
1388  $ rx(i,7,e)*zm(i,1)+rx(i,8,e)*zm(i,2)+rx(i,9,e)*zm(i,3)
1389  enddo
1390 
1391 
1392  call grad_rst(zr,zs,zt,zp(1,1),lxd,if3d)
1393  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1394  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1395  enddo
1396  call intp_rstd(fx,wk,lx1,lxd,if3d,1) ! Project back to coarse
1397 
1398  call grad_rst(zr,zs,zt,zp(1,2),lxd,if3d)
1399  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1400  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1401  enddo
1402  call intp_rstd(fy,wk,lx1,lxd,if3d,1) ! Project back to coarse
1403 
1404  call grad_rst(zr,zs,zt,zp(1,3),lxd,if3d)
1405  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1406  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1407  enddo
1408  call intp_rstd(fz,wk,lx1,lxd,if3d,1) ! Project back to coarse
1409 
1410 
1411  do i=1,nxyzd ! Convert convector (zp) to r-s-t coordinates
1412  tr(i,1)=
1413  $ rx(i,1,e)*zp(i,1)+rx(i,2,e)*zp(i,2)+rx(i,3,e)*zp(i,3)
1414  tr(i,2)=
1415  $ rx(i,4,e)*zp(i,1)+rx(i,5,e)*zp(i,2)+rx(i,6,e)*zp(i,3)
1416  tr(i,3)=
1417  $ rx(i,7,e)*zp(i,1)+rx(i,8,e)*zp(i,2)+rx(i,9,e)*zp(i,3)
1418  enddo
1419 
1420 
1421  call grad_rst(zr,zs,zt,zm(1,1),lxd,if3d)
1422  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1423  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1424  enddo
1425  call intp_rstd(gx,wk,lx1,lxd,if3d,1) ! Project back to coarse
1426 
1427  call grad_rst(zr,zs,zt,zm(1,2),lxd,if3d)
1428  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1429  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1430  enddo
1431  call intp_rstd(gy,wk,lx1,lxd,if3d,1) ! Project back to coarse
1432 
1433  call grad_rst(zr,zs,zt,zm(1,3),lxd,if3d)
1434  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1435  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)+tr(i,3)*zt(i)
1436  enddo
1437  call intp_rstd(gz,wk,lx1,lxd,if3d,1) ! Project back to coarse
1438 
1439  tmp = -0.5 ! vtrans() assumed to be 1.0 !
1440  do i=1,nxyz1
1441  bfx(i,1,1,e) = bfx(i,1,1,e) + tmp*( fx(i) + gx(i) )
1442  bmx(i,1,1,e) = bmx(i,1,1,e) + tmp*( fx(i) - gx(i) )
1443 
1444  bfy(i,1,1,e) = bfy(i,1,1,e) + tmp*( fy(i) + gy(i) )
1445  bmy(i,1,1,e) = bmy(i,1,1,e) + tmp*( fy(i) - gy(i) )
1446 
1447  bfz(i,1,1,e) = bfz(i,1,1,e) + tmp*( fz(i) + gz(i) )
1448  bmz(i,1,1,e) = bmz(i,1,1,e) + tmp*( fz(i) - gz(i) )
1449  enddo
1450 
1451  enddo
1452 
1453  else ! 2D
1454 c
1455  do e=1,nelv
1456 
1457 c Interpolate z+ and z- into fine mesh, translate to r-s-t coords
1458 
1459  call add3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
1460  call intp_rstd(zp(1,1),wk,lx1,lxd,if3d,0) ! 0 --> forward
1461 
1462  call add3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
1463  call intp_rstd(zp(1,2),wk,lx1,lxd,if3d,0)
1464 
1465  call sub3(wk,vx(1,1,1,e),bx(1,1,1,e),nxyz1)
1466  call intp_rstd(zm(1,1),wk,lx1,lxd,if3d,0)
1467 
1468  call sub3(wk,vy(1,1,1,e),by(1,1,1,e),nxyz1)
1469  call intp_rstd(zm(1,2),wk,lx1,lxd,if3d,0)
1470 
1471  do i=1,nxyzd ! Convert convector (zm) to r-s-t coordinates
1472  tr(i,1)=
1473  $ rx(i,1,e)*zm(i,1)+rx(i,2,e)*zm(i,2)
1474  tr(i,2)=
1475  $ rx(i,3,e)*zm(i,1)+rx(i,4,e)*zm(i,2)
1476  enddo
1477 
1478  call grad_rst(zr,zs,zt,zp(1,1),lxd,if3d)
1479  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1480  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
1481  enddo
1482  call intp_rstd(fx,wk,lx1,lxd,if3d,1) ! Project back to coarse
1483 
1484  call grad_rst(zr,zs,zt,zp(1,2),lxd,if3d)
1485  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1486  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
1487  enddo
1488  call intp_rstd(fy,wk,lx1,lxd,if3d,1) ! Project back to coarse
1489 
1490  do i=1,nxyzd ! Convert convector (zp) to r-s-t coordinates
1491  tr(i,1)=
1492  $ rx(i,1,e)*zp(i,1)+rx(i,2,e)*zp(i,2)
1493  tr(i,2)=
1494  $ rx(i,3,e)*zp(i,1)+rx(i,4,e)*zp(i,2)
1495  enddo
1496 
1497  call grad_rst(zr,zs,zt,zm(1,1),lxd,if3d)
1498  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1499  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
1500  enddo
1501  call intp_rstd(gx,wk,lx1,lxd,if3d,1) ! Project back to coarse
1502 
1503  call grad_rst(zr,zs,zt,zm(1,2),lxd,if3d)
1504  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
1505  wk(i) = tr(i,1)*zr(i)+tr(i,2)*zs(i)
1506  enddo
1507  call intp_rstd(gy,wk,lx1,lxd,if3d,1) ! Project back to coarse
1508 
1509  tmp = -0.5 ! vtrans() assumed to be 1.0 !
1510  do i=1,nxyz1
1511  bfx(i,1,1,e) = bfx(i,1,1,e) + tmp*( fx(i) + gx(i) )
1512  bmx(i,1,1,e) = bmx(i,1,1,e) + tmp*( fx(i) - gx(i) )
1513 
1514  bfy(i,1,1,e) = bfy(i,1,1,e) + tmp*( fy(i) + gy(i) )
1515  bmy(i,1,1,e) = bmy(i,1,1,e) + tmp*( fy(i) - gy(i) )
1516  enddo
1517  enddo
1518  endif
1519 
1520  return
1521  end
1522 c-----------------------------------------------------------------------
1523  subroutine cfl_check
1524 c
1525  include 'SIZE'
1526  include 'INPUT'
1527  include 'SOLN'
1528  include 'MASS'
1529  include 'TSTEP'
1530 c
1531  common /scrns/ ta1(lx1,ly1,lz1,lelv)
1532  $ , ta2(lx1,ly1,lz1,lelv)
1533  $ , ta3(lx1,ly1,lz1,lelv)
1534  $ , tb1(lx1,ly1,lz1,lelv)
1535  $ , tb2(lx1,ly1,lz1,lelv)
1536  $ , tb3(lx1,ly1,lz1,lelv)
1537 
1538 
1539  call opcopy(ta1,ta2,ta3,vx,vy,vz)
1540  call opcopy(tb1,tb2,tb3,bx,by,bz)
1541 
1542  ntot1 = lx1*ly1*lz1*nelv
1543  call phys_to_elsasser(ta1,ta2,ta3,tb1,tb2,tb3,ntot1) ! crude, but effective
1544 
1545  call compute_cfl(cflp,ta1,ta2,ta3,dt) ! vx = U+B
1546  call compute_cfl(cflm,tb1,tb2,tb3,dt) ! bx = U-B
1547 
1548  courno = max(cflp,cflm)
1549 
1550  if (nio.eq.0) write(6,1) istep,time,dt,cflp,cflm
1551  1 format(i9,1p4e15.7,' CFL')
1552 
1553  return
1554  end
1555 c--------------------------------------------------------------------
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
bcneutr
subroutine bcneutr
Definition: bdry.f:1200
exitt
void exitt()
Definition: comm_mpi.f:604
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
makeuq
subroutine makeuq
Definition: conduct.f:106
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
set_dealias_rx
subroutine set_dealias_rx
Definition: convect2.f:121
intp_rstd
subroutine intp_rstd(ju, u, mx, md, if3d, idir)
Definition: convect.f:347
grad_rst
subroutine grad_rst(ur, us, ut, u, md, if3d)
Definition: convect.f:615
vol_flow
subroutine vol_flow
Definition: drive2.f:1322
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
lagbfield
subroutine lagbfield
Definition: induct.f:55
makebsource_mhd
subroutine makebsource_mhd
Definition: induct.f:74
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
makebdfb
subroutine makebdfb
Definition: induct.f:179
phys_to_elsasser2
subroutine phys_to_elsasser2(u1, b1, n)
Definition: induct.f:824
specx
subroutine specx(b, nb, a, na, ba, ab, w)
Definition: induct.f:755
spec_curl_e
subroutine spec_curl_e(cb, b1, b2, b3, rx, ry, rz, sx, sy, sz, tx, ty, tz)
Definition: induct.f:687
advab_elsasser_fast
subroutine advab_elsasser_fast
Definition: induct.f:1326
elsasser_to_phys2
subroutine elsasser_to_phys2(u1, b1, n)
Definition: induct.f:838
cfl_check
subroutine cfl_check
Definition: induct.f:1524
extrapp
subroutine extrapp(p, plag)
Definition: induct.f:370
opzero
subroutine opzero(ux, uy, uz)
Definition: induct.f:406
opnorm
subroutine opnorm(unorm, ux, uy, uz, type3)
Definition: induct.f:419
elsasser_to_phys
subroutine elsasser_to_phys(u1, u2, u3, b1, b2, b3, n)
Definition: induct.f:797
setrhsp
subroutine setrhsp(p, h1, h2, h2inv, pset, niprev)
Definition: induct.f:1152
set_ifbcor
subroutine set_ifbcor
Definition: induct.f:1092
gensolnp
subroutine gensolnp(p, h1, h2, h2inv, pset, nprev)
Definition: induct.f:1222
econjp
subroutine econjp(pset, nprev, h1, h2, h2inv, ierr)
Definition: induct.f:1267
induct
subroutine induct(igeom)
Definition: induct.f:12
makextb
subroutine makextb
Definition: induct.f:140
lorentz_force
subroutine lorentz_force(lf, b1, b2, b3, w1, w2)
Definition: induct.f:447
cresvibp
subroutine cresvibp(resv1, resv2, resv3, h1, h2)
Definition: induct.f:252
cresvib
subroutine cresvib(resv1, resv2, resv3, h1, h2)
Definition: induct.f:222
compute_cfl
subroutine compute_cfl(cfl, u, v, w, dt)
Definition: induct.f:918
lorentz_force_e
subroutine lorentz_force_e(lf, b1, b2, b3, e)
Definition: induct.f:571
elsasserh
subroutine elsasserh(igeom)
Definition: induct.f:852
lorentz_force2
subroutine lorentz_force2(lf, b1, b2, b3)
Definition: induct.f:539
curl
subroutine curl(vort, u, v, w, ifavg, work1, work2)
Definition: induct.f:486
getdr
subroutine getdr(dri, zgm1, lx1)
Definition: induct.f:1009
incomprn
subroutine incomprn(ux, uy, uz, up)
Definition: induct.f:282
makeufb
subroutine makeufb
Definition: induct.f:123
phys_to_elsasser
subroutine phys_to_elsasser(u1, u2, u3, b1, b2, b3, n)
Definition: induct.f:770
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
add2col2
subroutine add2col2(a, b, c, n)
Definition: math.f:1565
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
add3
subroutine add3(a, b, c, n)
Definition: math.f:644
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
gllog
subroutine gllog(la, lb)
Definition: math.f:986
glsum
function glsum(x, n)
Definition: math.f:861
add3s2
subroutine add3s2(a, b, c, c1, c2, n)
Definition: math.f:716
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
invcol1
subroutine invcol1(a, n)
Definition: math.f:62
add2s1
subroutine add2s1(a, b, c1, n)
Definition: math.f:678
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
esolver
subroutine esolver(RES, H1, H2, H2INV, INTYPE)
Definition: navier0.f:2
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
makeuf
subroutine makeuf
Definition: navier1.f:1464
makebdf
subroutine makebdf
Definition: navier1.f:1549
specmp
subroutine specmp(b, nb, a, na, ba, ab, w)
Definition: navier1.f:3020
opcolv3c
subroutine opcolv3c(a1, a2, a3, b1, b2, b3, c, d)
Definition: navier1.f:2788
ortho
subroutine ortho(respr)
Definition: navier1.f:224
opbinv
subroutine opbinv(out1, out2, out3, inp1, inp2, inp3, h2inv)
Definition: navier1.f:776
opadd2col
subroutine opadd2col(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2751
opcopy
subroutine opcopy(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2474
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
opcolv2
subroutine opcolv2(a1, a2, a3, b1, b2)
Definition: navier1.f:2712
lagvel
subroutine lagvel
Definition: navier1.f:1930
setmap
subroutine setmap(n1, nd)
Definition: navier1.f:3051
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
setproj
subroutine setproj(n1, nd)
Definition: navier1.f:4117
makeabf
subroutine makeabf
Definition: navier1.f:1598
nekuf
subroutine nekuf(f1, f2, f3)
Definition: navier1.f:1483
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
dudxyz
subroutine dudxyz(du, u, rm1, sm1, tm1, jm1, imsh, isd)
Definition: navier1.f:1157
opsub2
subroutine opsub2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2360
hsolve
subroutine hsolve(name, u, r, h1, h2, vmk, vml, imsh, tol, maxit, isd, approx, napprox, bi)
Definition: navier4.f:564
chkptol
subroutine chkptol
Definition: navier4.f:269
vlsc3
function vlsc3(X, Y, B, N)
Definition: navier4.f:324
updrhse
subroutine updrhse(p, h1, h2, h2inv, ierr)
Definition: navier4.f:357
q_filter
subroutine q_filter(wght)
Definition: navier5.f:3
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
opadd2cm
subroutine opadd2cm(a1, a2, a3, b1, b2, b3, c)
Definition: plan4.f:439
unorm
subroutine unorm
Definition: subs1.f:352
hmhzsf
subroutine hmhzsf(name, u1, u2, u3, r1, r2, r3, h1, h2, rmask1, rmask2, rmask3, rmult, tol, maxit, matmod)
Definition: subs1.f:1096
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341