KTH framework for Nek5000 toolboxes; testing version  0.0.1
perturb.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine fluidp (igeom)
3 c
4 c Driver for perturbation velocity
5 c
6  include 'SIZE'
7  include 'INPUT'
8  include 'TSTEP'
9  include 'SOLN'
10 
11  do jp=1,npert
12 
13  if (nio.eq.0.and.igeom.eq.2) write(6,1) istep,time,jp
14  1 format(i9,1pe14.7,' Perturbation Solve:',i5)
15 
16  call perturbv (igeom)
17 
18  enddo
19 
20  jp=0 ! set jp to zero, for baseline flow
21 
22  return
23  end
24 c-----------------------------------------------------------------------
25  subroutine perturbv (igeom)
26 c
27 c
28 c Solve the convection-diffusion equation for the perturbation field,
29 c with projection onto a div-free space.
30 c
31 c
32  include 'SIZE'
33  include 'INPUT'
34  include 'EIGEN'
35  include 'SOLN'
36  include 'TSTEP'
37  include 'MASS'
38 C
39  COMMON /scrns/ resv1(lx1,ly1,lz1,lelv)
40  $ , resv2(lx1,ly1,lz1,lelv)
41  $ , resv3(lx1,ly1,lz1,lelv)
42  $ , dv1(lx1,ly1,lz1,lelv)
43  $ , dv2(lx1,ly1,lz1,lelv)
44  $ , dv3(lx1,ly1,lz1,lelv)
45  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
46  $ , h2(lx1,ly1,lz1,lelv)
47 c
48  ifield = 1
49 c
50  if (igeom.eq.1) then
51 c
52 c Old geometry, old velocity
53 c
54  call makefp
55  call lagfieldp
56 c
57  else
58 c
59 c New geometry, new velocity
60 c
61  intype = -1
62  call sethlm (h1,h2,intype)
63  call cresvipp (resv1,resv2,resv3,h1,h2)
64  call ophinv (dv1,dv2,dv3,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
65  call opadd2 (vxp(1,jp),vyp(1,jp),vzp(1,jp),dv1,dv2,dv3)
66  call incomprp (vxp(1,jp),vyp(1,jp),vzp(1,jp),prp(1,jp))
67 c
68  endif
69 c
70  return
71  end
72 c-----------------------------------------------------------------------
73  subroutine lagfieldp
74 c
75 c Keep old Vp-field(s)
76 c
77  include 'SIZE'
78  include 'INPUT'
79  include 'SOLN'
80  include 'TSTEP'
81 c
82  do ilag=nbdinp-1,2,-1
83  call opcopy
84  $ (vxlagp(1,ilag ,jp),vylagp(1,ilag ,jp),vzlagp(1,ilag ,jp)
85  $ ,vxlagp(1,ilag-1,jp),vylagp(1,ilag-1,jp),vzlagp(1,ilag-1,jp))
86  enddo
87  call opcopy(vxlagp(1,1,jp),vylagp(1,1,jp),vzlagp(1,1,jp)
88  $ ,vxp(1,jp) ,vyp(1,jp) ,vzp(1,jp) )
89 c
90  return
91  end
92 c-----------------------------------------------------------------------
93  subroutine makefp
94 c
95 c Make rhs for velocity perturbation equation
96 c
97  include 'SIZE'
98  include 'SOLN'
99  include 'MASS'
100  include 'INPUT'
101  include 'TSTEP'
102  include 'ADJOINT'
103 
104  call makeufp
105  if (filtertype.eq.2) call make_hpf
106  if (ifnav.and.(.not.ifchar).and.(.not.ifadj))call advabp
107  if (ifnav.and.(.not.ifchar).and.( ifadj))call advabp_adjoint
108  if (iftran) call makextp
109  call makebdfp
110 c
111  return
112  end
113 c--------------------------------------------------------------------
114  subroutine makeufp
115 c
116 c Compute and add: (1) user specified forcing function (FX,FY,FZ)
117 c
118  include 'SIZE'
119  include 'SOLN'
120  include 'MASS'
121  include 'TSTEP'
122 C
123  time = time-dt
124  call nekuf (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp))
125  call opcolv2 (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp)
126  $ ,vtrans(1,1,1,1,ifield),bm1)
127  time = time+dt
128 c
129  return
130  end
131 c--------------------------------------------------------------------
132  subroutine advabp
133 C
134 C Eulerian scheme, add convection term to forcing function
135 C at current time step.
136 C
137  include 'SIZE'
138  include 'INPUT'
139  include 'SOLN'
140  include 'MASS'
141  include 'TSTEP'
142 C
143  COMMON /scrns/ ta1(lx1*ly1*lz1*lelv)
144  $ , ta2(lx1*ly1*lz1*lelv)
145  $ , ta3(lx1*ly1*lz1*lelv)
146  $ , tb1(lx1*ly1*lz1*lelv)
147  $ , tb2(lx1*ly1*lz1*lelv)
148  $ , tb3(lx1*ly1*lz1*lelv)
149 C
150  ntot1 = lx1*ly1*lz1*nelv
151  ntot2 = lx2*ly2*lz2*nelv
152 c
153  if (if3d) then
154  call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
155  call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
156  call convop (ta1,tb1) ! du.grad U
157  call convop (ta2,tb2)
158  call convop (ta3,tb3)
159  call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
160 c
161  do i=1,ntot1
162  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
163  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
164  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
165  bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
166  enddo
167 c
168  call convop (ta1,vxp(1,jp)) ! U.grad dU
169  call convop (ta2,vyp(1,jp))
170  call convop (ta3,vzp(1,jp))
171 c
172  do i=1,ntot1
173  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
174  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
175  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
176  bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
177  enddo
178 c
179  else ! 2D
180 c
181  call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
182  call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
183  call convop (ta1,tb1) ! du.grad U
184  call convop (ta2,tb2)
185  call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
186 c
187  do i=1,ntot1
188  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
189  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
190  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
191  enddo
192 c
193  call convop (ta1,vxp(1,jp)) ! U.grad dU
194  call convop (ta2,vyp(1,jp))
195 c
196  do i=1,ntot1
197  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
198  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
199  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
200  enddo
201 c
202  endif
203 c
204  return
205  end
206 c--------------------------------------------------------------------
207  subroutine advabp_adjoint
208 C
209 C Eulerian scheme, add convection term to forcing function
210 C at current time step for backward part of adjoint:
211 C Convective term is now (U.Grad)u - (Grad U)^T .u
212 C instead of (u.Grad)U + (U.Grad)u in above subroutine advabp
213 C
214  include 'SIZE'
215  include 'INPUT'
216  include 'SOLN'
217  include 'MASS'
218  include 'TSTEP'
219  include 'GEOM'
220  include 'ADJOINT'
221 C
222  COMMON /scrns/ ta1(lx1*ly1*lz1*lelv)
223  $ , ta2(lx1*ly1*lz1*lelv)
224  $ , ta3(lx1*ly1*lz1*lelv)
225  $ , tb1(lx1*ly1*lz1*lelv)
226  $ , tb2(lx1*ly1*lz1*lelv)
227  $ , tb3(lx1*ly1*lz1*lelv)
228 
229 
230  real fact,factx,facty
231 C
232  ntot1 = lx1*ly1*lz1*nelv
233  ntot2 = lx2*ly2*lz2*nelv !dimensionn arrays
234  ntot = lx1*ly1*lz1*nelt
235 
236 c
237  if (if3d) then
238  call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
239  call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- u
240 c
241  call convop_adj (ta1,ta2,ta3,tb1,tb2,tb3,vx,vy,vz) ! u.grad U^T
242 
243  call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
244 c
245  do i=1,ntot1
246  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
247  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
248  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
249  bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
250  enddo
251 c
252 c
253 c
254  call convop (ta1,vxp(1,jp)) ! U.grad u
255  call convop (ta2,vyp(1,jp))
256  call convop (ta3,vzp(1,jp))
257 c
258  do i=1,ntot1
259  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
260  bfxp(i,jp) = bfxp(i,jp)+tmp*ta1(i)
261  bfyp(i,jp) = bfyp(i,jp)+tmp*ta2(i)
262  bfzp(i,jp) = bfzp(i,jp)+tmp*ta3(i)
263  enddo
264 c
265  if (ifheat) then ! dt.(grad T)
266 c ! term coming from the temperature convection
267  call opcolv3 (ta1,ta2,ta3,dtdx,dtdy,dtdz,tp)
268 c
269  do i=1,ntot1
270  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,2)
271  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
272  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
273  bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
274  enddo
275  endif
276 c
277  else ! 2D
278 
279  call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save velocity
280  call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
281 
282  call convop_adj (ta1,ta2,ta3,tb1,tb2,tb3,vx,vy,vz) ! u.((grad U)^T)
283 
284  call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity
285 
286  do i=1,ntot1
287  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
288  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
289  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
290  enddo
291 
292  call convop (ta1,vxp(1,jp)) ! U.grad u
293  call convop (ta2,vyp(1,jp))
294 c
295  do i=1,ntot1
296  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
297  bfxp(i,jp) = bfxp(i,jp)+tmp*ta1(i)
298  bfyp(i,jp) = bfyp(i,jp)+tmp*ta2(i)
299  enddo
300 
301  if (ifheat) then ! dt.(grad T)^T
302  ! term coming from the temperature convection
303  call opcolv3 (ta1,ta2,ta3,dtdx,dtdy,dtdz,tp)
304 c
305  do i=1,ntot1
306  tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,2)
307  bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
308  bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
309  enddo
310  endif
311 
312  endif
313 c
314  return
315  end
316 c--------------------------------------------------------------------
317  subroutine convop_adj (bdux,bduy,bduz,udx,udy,udz,cx,cy,cz)
318 
319  include 'SIZE'
320  include 'TOTAL'
321 
322  parameter(lxy=lx1*ly1*lz1,ltd=lxd*lyd*lzd)
323  common /scrcv/ fx(ltd),fy(ltd),fz(ltd)
324  $ , uf1(ltd),uf2(ltd),uf3(ltd),uf4(ltd),uf5(ltd),uf6(ltd)
325  real urx(ltd),usx(ltd),utx(ltd)
326  real ury(ltd),usy(ltd),uty(ltd)
327  real urz(ltd),usz(ltd),utz(ltd)
328  real bdux(1),bduy(1),bduz(1),u(1),cx(1),cy(1),cz(1)
329  real udx(1),udy(1),udz(1)
330  logical ifuf,ifcf ! u and/or c already on fine mesh?
331  integer e
332  real bdrx(1), bdry(1),bdrz (1)
333 
334  call set_dealias_rx
335  nxyz1 = lx1*ly1*lz1
336 c AM DING DING
337  nxyzd = lxd*lyd*lzd
338  nxyzu = nxyz1
339  nxyzc = nxyz1
340  ntot1=lx1*ly1*lz1*nelv
341  ic = 1 ! pointer to vector field C
342  iu = 1 ! pointer to scalar field u
343  ib = 1 ! pointer to scalar field Bdu
344  if(if3d)then
345  do e=1,nelv
346  ! map coarse velocity to fine mesh (C-->F)
347  call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
348  call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0)
349  call intp_rstd(fz,cz(ic),lx1,lxd,if3d,0)
350 
351  call intp_rstd(uf1,udx(iu),lx1,lxd,if3d,0) ! 0 --> forward
352  call grad_rst(urx,usx,utx,uf1,lxd,if3d)
353 
354  call intp_rstd(uf2,udy(iu),lx1,lxd,if3d,0)
355  call grad_rst(ury,usy,uty,uf2,lxd,if3d)
356 
357  call intp_rstd(uf3,udz(iu),lx1,lxd,if3d,0)
358  call grad_rst(urz,usz,utz,uf3,lxd,if3d)
359 
360  do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
361  uf4(i)=fx(i)*(rx(i,1,e)*urx(i)+rx(i,4,e)*usx(i)
362  $ +rx(i,7,e)*utx(i))+
363  $ fy(i)*(rx(i,1,e)*ury(i)+rx(i,4,e)*usy(i)
364  $ +rx(i,7,e)*uty(i))+
365  $ fz(i)*(rx(i,1,e)*urz(i)+rx(i,4,e)*usz(i)
366  $ +rx(i,7,e)*utz(i))
367  uf5(i)=fx(i)*(rx(i,2,e)*urx(i)+rx(i,5,e)*usx(i)
368  $ +rx(i,8,e)*utx(i))+
369  $ fy(i)*(rx(i,2,e)*ury(i)+rx(i,5,e)*usy(i)
370  $ +rx(i,8,e)*uty(i))+
371  $ fz(i)*(rx(i,2,e)*urz(i)+rx(i,5,e)*usz(i)
372  $ +rx(i,8,e)*utz(i))
373  uf6(i)=fx(i)*(rx(i,3,e)*urx(i)+rx(i,6,e)*usx(i)
374  $ +rx(i,9,e)*utx(i))+
375  $ fy(i)*(rx(i,3,e)*ury(i)+rx(i,6,e)*usy(i)
376  $ +rx(i,9,e)*uty(i))+
377  $ fz(i)*(rx(i,3,e)*urz(i)+rx(i,6,e)*usz(i)
378  $ +rx(i,9,e)*utz(i))
379  enddo
380 
381  call intp_rstd(bdux(ib),uf4,lx1,lxd,if3d,1) ! Project back to coarse
382  call intp_rstd(bduy(ib),uf5,lx1,lxd,if3d,1)
383  call intp_rstd(bduz(ib),uf6,lx1,lxd,if3d,1)
384 
385  ic = ic + nxyzc
386  iu = iu + nxyzu
387  ib = ib + nxyz1
388  enddo
389  call invcol2 (bdux,bm1,ntot1) ! local mass inverse
390  call invcol2 (bduy,bm1,ntot1) ! local mass inverse
391  call invcol2 (bduz,bm1,ntot1) ! local mass inverse
392  else
393  do e=1,nelv
394 
395  ! map coarse velocity to fine mesh (C-->F)
396  call intp_rstd(fx,cx(ic),lx1,lxd,if3d,0) ! 0 --> forward
397  call intp_rstd(fy,cy(ic),lx1,lxd,if3d,0)
398 
399  call intp_rstd(uf1,udx(iu),lx1,lxd,if3d,0)
400  call grad_rst(urx,usx,utx,uf1,lxd,if3d)
401 
402  call intp_rstd(uf2,udy(iu),lx1,lxd,if3d,0)
403  call grad_rst(ury,usy,uty,uf2,lxd,if3d)
404 
405  do i=1,nxyzd
406  uf4(i) = fx(i)*(rx(i,1,e)*urx(i)+rx(i,3,e)*usx(i))+
407  $ fy(i)*(rx(i,1,e)*ury(i)+rx(i,3,e)*usy(i))
408  uf5(i) = fx(i)*(rx(i,2,e)*urx(i)+rx(i,4,e)*usx(i))+
409  $ fy(i)*(rx(i,2,e)*ury(i)+rx(i,4,e)*usy(i))
410  enddo
411 
412  call intp_rstd(bdux(ib),uf4,lx1,lxd,if3d,1)
413  call intp_rstd(bduy(ib),uf5,lx1,lxd,if3d,1)
414 
415  ic = ic + nxyzc
416  iu = iu + nxyzu
417  ib = ib + nxyz1
418  end do
419  call invcol2 (bdux,bm1,ntot1) ! local mass inverse
420  call invcol2 (bduy,bm1,ntot1) ! local mass inverse
421  end if
422 
423 
424  return
425  end
426 c--------------------------------------------------------------------
427  subroutine makextp
428 c
429 c Add extrapolation terms to perturbation source terms
430 c
431 c (nek5 equivalent for velocity is "makeabf")
432 c
433  include 'SIZE'
434  include 'INPUT'
435  include 'SOLN'
436  include 'MASS'
437  include 'TSTEP'
438 C
439  common /scrns/ ta1(lx1,ly1,lz1,lelv)
440  $ , ta2(lx1,ly1,lz1,lelv)
441  $ , ta3(lx1,ly1,lz1,lelv)
442 c
443  ntot1 = lx1*ly1*lz1*nelv
444 c
445  ab0 = ab(1)
446  ab1 = ab(2)
447  ab2 = ab(3)
448  call add3s2 (ta1,exx1p(1,jp),exx2p(1,jp),ab1,ab2,ntot1)
449  call add3s2 (ta2,exy1p(1,jp),exy2p(1,jp),ab1,ab2,ntot1)
450  call copy (exx2p(1,jp),exx1p(1,jp),ntot1)
451  call copy (exy2p(1,jp),exy1p(1,jp),ntot1)
452  call copy (exx1p(1,jp),bfxp(1,jp),ntot1)
453  call copy (exy1p(1,jp),bfyp(1,jp),ntot1)
454  call add2s1 (bfxp(1,jp),ta1,ab0,ntot1)
455  call add2s1 (bfyp(1,jp),ta2,ab0,ntot1)
456  if (if3d) then
457  call add3s2 (ta3,exz1p(1,jp),exz2p(1,jp),ab1,ab2,ntot1)
458  call copy (exz2p(1,jp),exz1p(1,jp),ntot1)
459  call copy (exz1p(1,jp),bfzp(1,jp),ntot1)
460  call add2s1 (bfzp(1,jp),ta3,ab0,ntot1)
461  endif
462 c
463  return
464  end
465 c--------------------------------------------------------------------
466  subroutine makebdfp
467 C
468 C Add contributions to perturbation source from lagged BD terms.
469 C
470  include 'SIZE'
471  include 'SOLN'
472  include 'MASS'
473  include 'GEOM'
474  include 'INPUT'
475  include 'TSTEP'
476 C
477  COMMON /scrns/ ta1(lx1,ly1,lz1,lelv)
478  $ , ta2(lx1,ly1,lz1,lelv)
479  $ , ta3(lx1,ly1,lz1,lelv)
480  $ , tb1(lx1,ly1,lz1,lelv)
481  $ , tb2(lx1,ly1,lz1,lelv)
482  $ , tb3(lx1,ly1,lz1,lelv)
483  $ , h2(lx1,ly1,lz1,lelv)
484 C
485  ntot1 = lx1*ly1*lz1*nelv
486  const = 1./dt
487  call cmult2(h2,vtrans(1,1,1,1,ifield),const,ntot1)
488  call opcolv3c (tb1,tb2,tb3
489  $ ,vxp(1,jp),vyp(1,jp),vzp(1,jp),bm1,bd(2))
490 C
491  do ilag=2,nbd
492  if (ifgeom) then
493  call opcolv3c(ta1,ta2,ta3,vxlagp(1,ilag-1,jp),
494  $ vylagp(1,ilag-1,jp),
495  $ vzlagp(1,ilag-1,jp),
496  $ bm1lag(1,1,1,1,ilag-1),bd(ilag+1))
497  else
498  call opcolv3c(ta1,ta2,ta3,vxlagp(1,ilag-1,jp),
499  $ vylagp(1,ilag-1,jp),
500  $ vzlagp(1,ilag-1,jp),
501  $ bm1 ,bd(ilag+1))
502  endif
503  call opadd2 (tb1,tb2,tb3,ta1,ta2,ta3)
504  enddo
505  call opadd2col (bfxp(1,jp),bfyp(1,jp),bfzp(1,jp),tb1,tb2,tb3,h2)
506 c
507  return
508  end
509 c-----------------------------------------------------------------------
510  subroutine cresvipp(resv1,resv2,resv3,h1,h2)
511 c
512 c Account for inhomogeneous Dirichlet boundary contributions
513 c in rhs of perturbation eqn.
514 c n
515 c Also, subtract off best estimate of grad p
516 c
517  include 'SIZE'
518  include 'TOTAL'
519  real resv1 (lx1,ly1,lz1,1)
520  real resv2 (lx1,ly1,lz1,1)
521  real resv3 (lx1,ly1,lz1,1)
522  real h1 (lx1,ly1,lz1,1)
523  real h2 (lx1,ly1,lz1,1)
524  common /scruz/ w1(lx1,ly1,lz1,lelv)
525  $ , w2(lx1,ly1,lz1,lelv)
526  $ , w3(lx1,ly1,lz1,lelv)
527  $ , prextr(lx2,ly2,lz2,lelv)
528 c
529  ntot1 = lx1*ly1*lz1*nelv
530  ntot2 = lx2*ly2*lz2*nelv
531 c
532  call bcdirvc (vxp(1,jp),vyp(1,jp),vzp(1,jp)
533  $ ,v1mask,v2mask,v3mask)
534  call extrapprp(prextr)
535  call opgradt(resv1,resv2,resv3,prextr)
536  call opadd2(resv1,resv2,resv3,bfxp(1,jp),bfyp(1,jp),bfzp(1,jp))
537  call ophx (w1,w2,w3,vxp(1,jp),vyp(1,jp),vzp(1,jp),h1,h2)
538  call opsub2(resv1,resv2,resv3,w1,w2,w3)
539 c
540  return
541  end
542 c--------------------------------------------------------------------
543  subroutine heatp (igeom)
544 C
545 C Driver for temperature or passive scalar.
546 C
547 C Current version:
548 C (1) Varaiable properties.
549 C (2) Implicit time stepping.
550 C (3) User specified tolerance for the Helmholtz solver
551 C (not based on eigenvalues).
552 C (4) A passive scalar can be defined on either the
553 C temperatur or the velocity mesh.
554 C (5) A passive scalar has its own multiplicity (B.C.).
555 C
556  include 'SIZE'
557  include 'INPUT'
558  include 'TSTEP'
559  include 'SOLN'
560 C
561  do jp=1,npert
562  do ifield=2,nfield
563  intype = -1
564  IF (.NOT.iftmsh(ifield)) imesh = 1
565  IF ( iftmsh(ifield)) imesh = 2
566  CALL unorm
567  CALL settolt
568  CALL cdscalp (igeom)
569  enddo
570  enddo
571 c
572  jp=0 ! set jp to zero, for baseline flow
573 c
574  return
575  end
576 C
577 C-----------------------------------------------------------------------
578  subroutine cdscalp (igeom)
579 C-----------------------------------------------------------------------
580 C
581 C Solve the convection-diffusion equation for passive scalar IPSCAL
582 C
583 C-----------------------------------------------------------------------
584  include 'SIZE'
585  include 'INPUT'
586  include 'GEOM'
587  include 'MVGEOM'
588  include 'SOLN'
589  include 'MASS'
590  include 'TSTEP'
591  COMMON /cprint/ ifprint
592  LOGICAL IFPRINT
593  LOGICAL IFCONV
594 C
595  COMMON /scrns/ ta(lx1,ly1,lz1,lelt)
596  $ ,tb(lx1,ly1,lz1,lelt)
597  COMMON /scrvh/ h1(lx1,ly1,lz1,lelt)
598  $ ,h2(lx1,ly1,lz1,lelt)
599 c
600  include 'ORTHOT'
601  ifld1 = ifield-1
602  napproxt(1,ifld1) = laxtt
603 C
604  IF (igeom.EQ.1) THEN
605 C
606 C Old geometry
607 C
608  CALL makeqp
609  CALL lagscalp
610 C
611  ELSE
612 C
613  IF (ifprint) THEN
614  IF (ifield.EQ.2.AND.nio.EQ.0) THEN
615  WRITE (6,*) ' Temperature/Passive scalar solution'
616  ENDIF
617  ENDIF
618  if1=ifield-1
619  write(name4t,1) if1
620  1 format('TEM',i1)
621 C
622 C New geometry
623 C
624  nel = nelfld(ifield)
625  ntot = lx1*ly1*lz1*nel
626 C
627  intype = 0
628  IF (iftran) intype = -1
629  CALL sethlm (h1,h2,intype)
630  CALL bcneusc (ta,-1)
631  CALL add2 (h2,ta,ntot)
632  CALL bcdirsc (tp(1,ifield-1,jp))
633  CALL axhelm (ta,tp(1,ifield-1,jp),h1,h2,imesh,1)
634  CALL sub3 (tb,bqp(1,ifield-1,jp),ta,ntot)
635  CALL bcneusc (ta,1)
636  CALL add2 (tb,ta,ntot)
637 c
638  CALL hmholtz (name4t,ta,tb,h1,h2
639  $ ,tmask(1,1,1,1,ifield-1)
640  $ ,tmult(1,1,1,1,ifield-1)
641  $ ,imesh,tolht(ifield),nmxt(ifield-1),1)
642 c
643  CALL add2 (tp(1,ifield-1,jp),ta,ntot)
644 C
645  CALL bcneusc (ta,1)
646  CALL add2 (bqp(1,ifield-1,jp),ta,ntot)
647 C
648  ENDIF
649 C
650  return
651  end
652 C
653  subroutine makeqp
654 C----------------------------------------------------------------------
655 C
656 C Generate forcing function for the solution of a passive scalar.
657 C !! NOTE: Do not change the content of the array BQ until the current
658 C time step is completed
659 C
660 C----------------------------------------------------------------------
661  include 'SIZE'
662  include 'GEOM'
663  include 'INPUT'
664  include 'TSTEP'
665  include 'SOLN'
666  CALL makeuqp
667  IF (filtertype.EQ.2) CALL make_hpf
668  IF (ifadvc(ifield).AND.(.NOT.ifchar)) CALL convabp
669  IF (iftran) CALL makeabqp
670  IF ((iftran.AND..NOT.ifchar).OR.
671  $ (iftran.AND..NOT.ifadvc(ifield).AND.ifchar)) CALL makebdqp
672 c IF (IFADVC(IFIELD).AND.IFCHAR) CALL CONVCHP
673  IF (ifadvc(ifield).AND.ifchar) write(6,*) 'no convchp'
674  IF (ifadvc(ifield).AND.ifchar) call exitt
675 c
676  return
677  end
678 C
679  subroutine makeuqp
680 C---------------------------------------------------------------------
681 C
682 C Fill up user defined forcing function and collocate will the
683 C mass matrix on the Gauss-Lobatto mesh.
684 C
685 C---------------------------------------------------------------------
686  include 'SIZE'
687  include 'MASS'
688  include 'SOLN'
689  include 'TSTEP'
690 C
691  ntot = lx1*ly1*lz1*nelfld(ifield)
692 C
693  time = time-dt ! time is tn
694 c
695  call rzero ( bqp(1,ifield-1,jp) , ntot)
696  call setqvol ( bqp(1,ifield-1,jp) )
697  call col2 ( bqp(1,ifield-1,jp) ,bm1,ntot)
698 c
699  time = time+dt ! restore time
700 C
701  return
702  end
703 C
704  subroutine convabp
705 C---------------------------------------------------------------
706 C
707 C Eulerian scheme, add convection term to forcing function
708 C at current time step.
709 C
710 C---------------------------------------------------------------
711  include 'SIZE'
712  include 'ADJOINT'
713  include 'SOLN'
714  include 'MASS'
715  include 'TSTEP'
716 c
717  common /scruz/ ta(lx1,ly1,lz1,lelt)
718  $ , ua(lx1,ly1,lz1,lelt)
719  $ , ub(lx1,ly1,lz1,lelt)
720  $ , uc(lx1,ly1,lz1,lelt)
721  real coeff
722 c
723  nel = nelfld(ifield)
724  ntot1 = lx1*ly1*lz1*nel
725 c
726  if (.not.ifadj) then
727  call opcopy(ua,ub,uc,vx,vy,vz)
728  call opcopy(vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
729  call convop(ta,t(1,1,1,1,ifield-1)) ! dU.grad T
730  call opcopy(vx,vy,vz,ua,ub,uc)
731  endif
732 c
733  call convop (ua,tp(1,ifield-1,jp)) ! U.grad dT
734 c
735  if (.not.ifadj) then
736  call add2 (ta,ua,ntot1) ! U.grad dT + dU.grad T
737  else
738  call copy (ta,ua,ntot1)
739  coeff=-1.0
740  call cmult(ta,coeff,ntot1) ! -U.grad dT
741  ! the second term depends on the buoyancy
742  endif
743 c
744  call col2 (ta,vtrans(1,1,1,1,ifield),ntot1) !vtrans (U.grad dT + dU.grad T)
745  call subcol3 (bqp(1,ifield-1,jp),bm1,ta,ntot1)
746 c
747  return
748  end
749 C
750  subroutine makeabqp
751 C-----------------------------------------------------------------------
752 C
753 C Sum up contributions to 3rd order Adams-Bashforth scheme.
754 C
755 C-----------------------------------------------------------------------
756  include 'SIZE'
757  include 'SOLN'
758  include 'TSTEP'
759 C
760  COMMON /scruz/ ta(lx1,ly1,lz1,lelt)
761 C
762  ab0 = ab(1)
763  ab1 = ab(2)
764  ab2 = ab(3)
765  nel = nelfld(ifield)
766  ntot1 = lx1*ly1*lz1*nel
767 C
768  CALL add3s2 (ta,vgradt1p(1,ifield-1,jp),
769  $ vgradt2p(1,ifield-1,jp),ab1,ab2,ntot1)
770  CALL copy ( vgradt2p(1,ifield-1,jp),
771  $ vgradt1p(1,ifield-1,jp),ntot1)
772  CALL copy ( vgradt1p(1,ifield-1,jp),
773  $ bqp(1,ifield-1,jp),ntot1)
774  CALL add2s1 (bqp(1,ifield-1,jp),ta,ab0,ntot1)
775 C
776  return
777  end
778 C
779  subroutine makebdqp
780 C-----------------------------------------------------------------------
781 C
782 C Add contributions to Q from lagged BD terms.
783 C
784 C-----------------------------------------------------------------------
785  include 'SIZE'
786  include 'SOLN'
787  include 'MASS'
788  include 'GEOM'
789  include 'INPUT'
790  include 'TSTEP'
791 C
792  COMMON /scrns/ ta(lx1,ly1,lz1,lelt)
793  $ , tb(lx1,ly1,lz1,lelt)
794  $ , h2(lx1,ly1,lz1,lelt)
795 C
796  nel = nelfld(ifield)
797  ntot1 = lx1*ly1*lz1*nel
798  const = 1./dt
799  CALL copy (h2,vtrans(1,1,1,1,ifield),ntot1)
800  CALL cmult (h2,const,ntot1)
801 C
802  CALL col3 (tb,bm1,tp(1,ifield-1,jp),ntot1)
803  CALL cmult (tb,bd(2),ntot1)
804 C
805  DO 100 ilag=2,nbd
806  IF (ifgeom) THEN
807  CALL col3 (ta,bm1lag(1,1,1,1,ilag-1),
808  $ tlagp(1,ilag-1,ifield-1,jp),ntot1)
809  ELSE
810  CALL col3 (ta,bm1,
811  $ tlagp(1,ilag-1,ifield-1,jp),ntot1)
812  ENDIF
813  CALL add2s2(tb,ta,bd(ilag+1),ntot1)
814  100 CONTINUE
815 C
816  CALL col2 (tb,h2,ntot1)
817  CALL add2 (bqp(1,ifield-1,jp),tb,ntot1)
818 C
819  return
820  end
821 C
822  subroutine lagscalp
823 C-----------------------------------------------------------------------
824 C
825 C Keep old passive scalar field(s)
826 C
827 C-----------------------------------------------------------------------
828  include 'SIZE'
829  include 'INPUT'
830  include 'SOLN'
831  include 'TSTEP'
832 C
833  ntot1 = lx1*ly1*lz1*nelfld(ifield)
834 C
835  DO 100 ilag=nbdinp-1,2,-1
836  CALL copy (tlagp(1,ilag ,ifield-1,jp),
837  $ tlagp(1,ilag-1,ifield-1,jp),ntot1)
838  100 CONTINUE
839 C
840  CALL copy (tlagp(1,1,ifield-1,jp),tp(1,ifield-1,jp),ntot1)
841 C
842  return
843  end
844 c-----------------------------------------------------------------------
845  subroutine incomprp (ux,uy,uz,up)
846 c
847 c Project U onto the closest incompressible field
848 c
849 c Input: U := (ux,uy,uz)
850 c
851 c Output: updated values of U, iproj, proj; and
852 c up := pressure correction req'd to impose div U = 0
853 c
854 c
855 c Dependencies: ifield ==> which "density" (vtrans) is used.
856 c
857 c
858  include 'SIZE'
859  include 'TOTAL'
860  include 'CTIMER'
861 c
862  common /scrns/ w1(lx1,ly1,lz1,lelv)
863  $ , w2(lx1,ly1,lz1,lelv)
864  $ , w3(lx1,ly1,lz1,lelv)
865  $ , dv1(lx1,ly1,lz1,lelv)
866  $ , dv2(lx1,ly1,lz1,lelv)
867  $ , dv3(lx1,ly1,lz1,lelv)
868  $ , dp(lx2,ly2,lz2,lelv)
869  common /scrvh/ h1(lx1,ly1,lz1,lelv)
870  $ , h2(lx1,ly1,lz1,lelv)
871  common /scrhi/ h2inv(lx1,ly1,lz1,lelv)
872  COMMON /scrch/ prextr(lx2,ly2,lz2,lelv)
873  logical ifprjp
874 
875 c
876  if (icalld.eq.0) tpres=0.0
877  icalld=icalld+1
878  npres=icalld
879 c
880  ntot1 = lx1*ly1*lz1*nelv
881  ntot2 = lx2*ly2*lz2*nelv
882  intype = 1
883  dtbd = bd(1)/dt
884 
885  call rzero (h1,ntot1)
886 c call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
887  call cmult2 (h2,vtrans(1,1,1,1,ifield),dtbd,ntot1)
888  call invers2 (h2inv,h2,ntot1)
889 
890  call opdiv (dp,ux,uy,uz)
891  call chsign (dp,ntot2)
892  call ortho (dp)
893 
894 
895 C******************************************************************
896 
897 
898  ifprjp=.false. ! project out previous pressure solutions?
899  istart=param(95)
900  if (istep.ge.istart.and.istart.ne.0) ifprjp=.true.
901 
902  ! Most likely, the following can be commented out. (pff, 1/6/2010)
903 c if (npert.gt.1.or.ifbase) ifprjp=.false.
904 cpff if (ifprjp) call setrhs (dp,h1,h2,h2inv)
905 
906  call esolver (dp,h1,h2,h2inv,intype)
907 
908 cpff if (ifprjp) call gensoln (dp,h1,h2,h2inv)
909 
910 cNOTE: The "cpff" comments added 11/24/17 to avoid old-style projection,
911 cNOTE: which should be replaced with something more updated.
912 
913 C******************************************************************
914 
915  call opgradt (w1 ,w2 ,w3 ,dp)
916  call opbinv (dv1,dv2,dv3,w1 ,w2 ,w3 ,h2inv)
917  call opadd2 (ux ,uy ,uz ,dv1,dv2,dv3)
918 
919  call extrapprp(prextr)
920  call lagpresp
921  call add3(up,prextr,dp,ntot2)
922 
923  return
924  end
925 c------------------------------------------------------------------------
926  subroutine extrapprp (prextr)
927 C
928 C Pressure extrapolation
929 C
930  include 'SIZE'
931  include 'SOLN'
932  include 'TSTEP'
933  COMMON /ctmp0/ dpr(lx2,ly2,lz2,lelv)
934  REAL PREXTR (LX2,LY2,LZ2,LELV)
935 C
936  ntot2 = lx2*ly2*lz2*nelv
937  if (nbd.eq.1.or.nbd.eq.2) then
938  call copy (prextr,prp(1,jp),ntot2)
939  elseif (nbd.eq.3) then
940  const = dtlag(1)/dtlag(2)
941  call sub3 (dpr,prp(1,jp),prlagp(1,1,jp),ntot2)
942  call cmult(dpr,const,ntot2)
943  call add3 (prextr,prp(1,jp),dpr,ntot2)
944  elseif (nbd.gt.3) then
945  write (6,*) 'Pressure extrapolation cannot be completed'
946  write (6,*) 'Try a lower-order temporal scheme'
947  call exitt
948  endif
949  return
950  end
951 C-------------------------------------------------------------------
952  subroutine lagpresp
953 C
954 C save old pressure values
955 C
956  include 'SIZE'
957  include 'SOLN'
958  include 'TSTEP'
959 C
960  if (nbdinp.eq.3) then
961  ntot2 = lx2*ly2*lz2*nelv
962  call copy (prlagp(1,1,jp),prp(1,jp),ntot2)
963  endif
964  return
965  end
966 C-------------------------------------------------------------------
967  subroutine lyap_scale ! Rescale / orthogonalize perturbation fields
968 c
969 c
970  include 'SIZE'
971  include 'TOTAL'
972 c
973  real sigma(0:lpert)
974 c
975  ntotv = lx1*ly1*lz1*nelv
976  ntotp = lx2*ly2*lz2*nelv
977  ntott = lx1*ly1*lz1*nelt
978 c
979  do j=1,npert
980  call normvc(h1,semi,pl2,plinf,vxp(1,j),vyp(1,j),vzp(1,j))
981  sigma(j) = pl2
982  if (pl2.gt.0) then
983  write(6,*) 'this is pl2:',pl2
984  scale = 1./pl2
985  call opcmult(vxp(1,j),vyp(1,j),vzp(1,j),scale)
986  call cmult(tp(1,1,j),scale,ntott)
987  call cmult(prp(1,j) ,scale,ntotp)
988  endif
989 c
990 c Have to do lag terms as well, etc
991 c
992 c Also, must orthogonalize
993  enddo
994 c
995  if (nio.eq.0) then
996  if (npert.eq.1) write(6,1) istep,time,(sigma(j),j=1,npert)
997  if (npert.eq.2) write(6,2) istep,time,(sigma(j),j=1,npert)
998  if (npert.eq.3) write(6,3) istep,time,(sigma(j),j=1,npert)
999  if (npert.eq.4) write(6,4) istep,time,(sigma(j),j=1,npert)
1000  if (npert.eq.5) write(6,5) istep,time,(sigma(j),j=1,npert)
1001  if (npert.eq.6) write(6,6) istep,time,(sigma(j),j=1,npert)
1002  if (npert.eq.7) write(6,7) istep,time,(sigma(j),j=1,npert)
1003  if (npert.eq.8) write(6,8) istep,time,(sigma(j),j=1,npert)
1004  if (npert.eq.9) write(6,9) istep,time,(sigma(j),j=1,npert)
1005  endif
1006  1 format(i8,1p2e12.4,' pgrow')
1007  2 format(i8,1p3e12.4,' pgrow')
1008  3 format(i8,1p4e12.4,' pgrow')
1009  4 format(i8,1p5e12.4,' pgrow')
1010  5 format(i8,1p6e12.4,' pgrow')
1011  6 format(i8,1p7e12.4,' pgrow')
1012  7 format(i8,1p8e12.4,' pgrow')
1013  8 format(i8,1p9e12.4,' pgrow')
1014  9 format(i8,1p10e12.4,' pgrow')
1015 c
1016  return
1017  end
1018 c-----------------------------------------------------------------------
1019  subroutine out_pert ! dump perturbation .fld files
1020 c
1021  include 'SIZE'
1022  include 'TOTAL'
1023  include 'NEKUSE'
1024 c
1025  character*3 s3
1026 c
1027  do jpp=1,npert
1028  write(s3,3) jpp
1029  3 format('p',i2.2)
1030  call outpost2
1031  $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),prp(1,jpp),tp(1,1,jpp),1,s3)
1032 c call writehist
1033 c $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),jpp)
1034  enddo
1035 c
1036  return
1037  end
1038 c-----------------------------------------------------------------------
1039  subroutine pert_add2s2(i,j,scale) ! xi = xi + scale * xj
1040  include 'SIZE'
1041  include 'TOTAL'
1042 
1043  ntotp = lx2*ly2*lz2*nelv
1044  ntotv = lx1*ly1*lz1*nelv
1045  ntott = lx1*ly1*lz1*nelt
1046 
1047  call add2s2(vxp(1,i),vxp(1,j),scale,ntotv)
1048  call add2s2(vyp(1,i),vyp(1,j),scale,ntotv)
1049  if (if3d) call add2s2(vzp(1,i),vzp(1,j),scale,ntotv)
1050  call add2s2(prp(1,i),prp(1,j),scale,ntotp)
1051 
1052  do l=1,lorder-1
1053  call add2s2(vxlagp(1,l,i),vxlagp(1,l,j),scale,ntotv)
1054  call add2s2(vylagp(1,l,i),vylagp(1,l,j),scale,ntotv)
1055  if (if3d) call add2s2(vzlagp(1,l,i),vzlagp(1,l,j),scale,ntotv)
1056  if (l.le.lorder-2)
1057  $ call add2s2(prlagp(1,l,i),prlagp(1,l,j),scale,ntotp)
1058  enddo
1059 
1060  call add2s2(exx1p(1,i),exx1p(1,j),scale,ntotv)
1061  call add2s2(exy1p(1,i),exy1p(1,j),scale,ntotv)
1062  if (if3d) call add2s2(exz1p(1,i),exz1p(1,j),scale,ntotv)
1063 
1064  call add2s2(exx2p(1,i),exx2p(1,j),scale,ntotv)
1065  call add2s2(exy2p(1,i),exy2p(1,j),scale,ntotv)
1066  if (if3d) call add2s2(exz2p(1,i),exz2p(1,j),scale,ntotv)
1067 
1068  if (ifheat) then
1069  do k=0,npscal
1070  k1=k+1
1071  ntotk = lx1*ly1*lz1*nelfld(k+2)
1072  call add2s2(tp(1,k1,i),tp(1,k1,j),scale,ntotk)
1073  do l=1,lorder-1
1074  call add2s2(tlagp(1,l,k1,i),tlagp(1,l,k1,j),scale,ntotk)
1075  enddo
1076  call add2s2(vgradt1p(1,k1,i),vgradt1p(1,k1,j),scale,ntotk)
1077  call add2s2(vgradt2p(1,k1,i),vgradt2p(1,k1,j),scale,ntotk)
1078  enddo
1079  endif
1080 
1081  return
1082  end
1083 c-----------------------------------------------------------------------
1084  function pert_inner_prod(i,j) ! weighted inner product vi^T vj
1085  include 'SIZE'
1086  include 'TOTAL'
1087 
1088  common/normset/pran, ray, rayc
1089 
1090  ntotv=lx1*ly1*lz1*nelv
1091  ntott=lx1*ly1*lz1*nelt
1092 
1093  s1 = rayc*glsc3(vxp(1,i),bm1,vxp(1,j),ntotv)
1094  s2 = rayc*glsc3(vyp(1,i),bm1,vyp(1,j),ntotv)
1095  s3 = 0
1096  if (if3d) s3 = rayc*glsc3(vzp(1,i),bm1,vzp(1,j),ntotv)
1097 
1098  t1 = 0
1099  if (ifheat) t1=pran*ray*ray*glsc3(tp(1,1,i),bm1,tp(1,1,j),ntott)
1100 
1101  pert_inner_prod = (s1+s2+s3+t1)/volvm1
1102 
1103  return
1104  end
1105 c-----------------------------------------------------------------------
1106  subroutine pert_ortho_norm ! orthogonalize and rescale pert. arrays
1107  include 'SIZE'
1108  include 'TOTAL'
1109 
1110  do k=1,npert
1111  call pert_ortho_norm1(k)
1112  enddo
1113 
1114  return
1115  end
1116 c-----------------------------------------------------------------------
1117  subroutine pert_ortho_norm1 (k) ! orthogonalize k against 1...k-1
1118  include 'SIZE'
1119  include 'TOTAL'
1120 
1121  do j=1,k-1
1122  scale = -pert_inner_prod(j,k)
1123  call pert_add2s2(k,j,scale) ! xi = xi + scale * xj
1124  enddo
1125  scale = pert_inner_prod(k,k)
1126  if (scale.gt.0) scale = 1./scale
1127  call rescalepert(pertnorm,scale,k)
1128 
1129  return
1130  end
1131 c-----------------------------------------------------------------------
1132  function opnorm2(v1,v2,v3)
1133  include 'SIZE'
1134  include 'TOTAL'
1135 c
1136  real v1(1) , v2(1), v3(1)
1137  real normsq1,normsq2,normsq3,opnorm
1138 c
1139  ntotv=lx1*ly1*lz1*nelv
1140  normsq1=glsc3(v1,bm1,v1,ntotv)
1141  normsq2=glsc3(v2,bm1,v2,ntotv)
1142  if(if3d) then
1143  normsq3=glsc3(v3,bm1,v3,ntotv)
1144  else
1145  normsq3=0
1146  endif
1147 
1148  opnorm2=normsq1+normsq2+normsq3
1149  if (opnorm2.gt.0) opnorm2=sqrt(opnorm2/volvm1)
1150  return
1151  end
1152 c-----------------------------------------------------------------------
1153 
1154  function tnorm(temp)
1155  include 'SIZE'
1156  include 'TOTAL'
1157 
1158  real temp(*)
1159 c
1160  ntotv = lx1*ly1*lz1*nelv
1161  tnorm = sqrt( glsc3(temp,bm1,temp,ntotv) /voltm1)
1162 c
1163  return
1164  end
1165 c--------------------------------------------
1166  function dmnorm1(v1,v2,v3,temp)
1167 c Norm weighted by mass matrix
1168  include 'SIZE'
1169  include 'TOTAL'
1170 
1171  real v1(1),v2(1),v3(1),temp(1)
1172  real normsq1,normsq2,normsq3,normsqt,dmnorm
1173  common/normset/pran, ray, rayc
1174 
1175  ntotv=lx1*ly1*lz1*nelv
1176  normsq1=(rayc)*glsc3(v1,bm1,v1,ntotv)
1177  normsq2=(rayc)*glsc3(v2,bm1,v2,ntotv)
1178  if(if3d) then
1179  normsq3=(rayc)*glsc3(v3,bm1,v3,ntotv)
1180  else
1181  normsq3=0
1182  endif
1183 
1184  if(ifheat) then
1185  normsqt = (pran*ray*ray)*glsc3(temp,bm1,temp,ntotv)
1186  else
1187  normsqt = 0
1188  endif
1189 
1190  dmnorm1=sqrt((normsq1+normsq2+normsq3+normsqt)/volvm1)
1191 
1192  return
1193  end
1194 
1195 c---------------------------------------------------------------
1196  subroutine opscale(v1,v2,v3,temp,alpha)
1197 c v = alpha*v
1198  include 'SIZE'
1199  include 'INPUT'
1200 
1201  real alpha
1202  real v1(1),v2(1),v3(1),temp(1)
1203 
1204  ltotv=lx1*ly1*lz1*lelv
1205  ltott=lx1*ly1*lz1*lelt
1206 
1207  call cmult(v1,alpha,ltotv)
1208  call cmult(v2,alpha,ltotv)
1209  if (if3d) call cmult(v3,alpha,ltotv)
1210  if (ifheat) call cmult(temp,alpha,ltott*ldimt)
1211 
1212  return
1213  end
1214 
1215 c---------------------------------------------------------------
1216  subroutine opscalev(v1,v2,v3,alpha)
1217 c v = alpha*v
1218  include 'SIZE'
1219  include 'INPUT'
1220  real alpha
1221  real v1(*),v2(*),v3(*)
1222 
1223  ntotv=lx1*ly1*lz1*nelv
1224 
1225  call cmult(v1,alpha,ntotv)
1226  call cmult(v2,alpha,ntotv)
1227 
1228  if (if3d) call cmult(v3,alpha,ntotv)
1229 c
1230  return
1231  end
1232 
1233 c-----------------------------------------------------------------------
1234  subroutine computelyap
1235  include 'SIZE'
1236  include 'TOTAL'
1237 
1238  do jpp=1,npert ! Loop through each Lyapunov eigenvalue
1239  call computelyap1
1240  $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),jpp)
1241  enddo
1242 
1243  return
1244  end
1245 c-----------------------------------------------------------------------
1246  subroutine computelyap1(vxq,vyq,vzq,tq,jpp)
1247  include 'SIZE'
1248  include 'TOTAL'
1249 
1250  real vxq(1),vyq(1),vzq(1),tq(1)
1251  real lyapsum,twt
1252  common /pertsave/ timestart,tinitial
1253 
1254  integer icount
1255  save icount
1256  data icount /0/
1257 
1258  logical if_restart,if_ortho_lyap
1259  common/restar/if_restart,if_ortho_lyap
1260 
1261  character*132 lyprestart
1262  common/restflename/lyprestart !file for restart data
1263 
1264  twt = param(126) !time to wait to start computing exponents
1265 
1266  if (nio.eq.0)
1267  $ write(6,8) istep,icount,time,twt,(lyap(k,jpp),k=1,3),jpp
1268  8 format(i9,i4,2f8.4,1p3e12.4,i3,'clyap')
1269 
1270  if(time.lt.twt) then
1271 c
1272 c For a fresh simulation, then we wait 5 vertical diffusion
1273 c times before we start measuring, so in this case just rescale
1274 c
1275  pertnorm = dmnorm1(vxq,vyq,vzq,tq)
1276  pertinvnorm = 1.0/pertnorm
1277  call rescalepert(pertnorm,pertinvnorm,jpp)
1278  lyap(3,jpp) = pertnorm !store latest norm
1279  timestart = time
1280  tinitial = time
1281  icount = 0
1282  return
1283  else
1284  if (jpp.eq.1) icount = icount + 1
1285  endif
1286 
1287  irescale = param(128)
1288  if (icount.eq.irescale) then
1289 
1290  lyapsum = lyap(2,jpp)
1291  oldpertnorm = lyap(3,jpp)
1292  pertnorm=dmnorm1(vxq,vyq,vzq,tq)
1293 c
1294  lyap(1,jpp) = log(pertnorm/oldpertnorm)/(time-timestart)
1295  lyapsum = lyapsum + log(pertnorm/oldpertnorm)
1296  lyap(2,jpp) = lyapsum
1297 
1298  if(nid.eq.0) then ! write out results to the .lyp file
1299 
1300  write(6 ,1) istep,time,lyap(1,jpp),lyapsum,pertnorm,jpp
1301  write(79,2) time,lyap(1,jpp),lyapsum,pertnorm,oldpertnorm,jpp
1302  1 format(i9,1p4e17.8,i4,'lyap')
1303  2 format(1p5e17.8,i4,'lyap')
1304 c call flushbuffer(79)
1305 
1306  if (jpp.eq.1) open(unit=96,file=lyprestart)
1307  write(96,9) lyapsum,timestart,jpp
1308  9 format(1p2e19.11,i9)
1309  if (jpp.eq.npert) close(96)
1310  endif
1311 
1312  pertinvnorm = 1.0/pertnorm
1313  call rescalepert(pertnorm,pertinvnorm,jpp)
1314  lyap(3,jpp) = pertnorm !save current pertnorm as old pertnorm
1315 
1316  if (jpp.eq.npert) then
1317  icount = 0
1318  timestart = time
1319  endif
1320 
1321  if_ortho_lyap = .false.
1322  if (param(125).gt.0) if_ortho_lyap = .true.
1323  if (jpp.eq.npert .and. if_ortho_lyap) call pert_ortho_norm
1324 
1325  endif
1326 
1327  return
1328  end
1329 c-----------------------------------------------------------------------
1330 
1331  subroutine rescalepert(pertnorm,pertinvnorm,jpp)
1332  include 'SIZE'
1333  include 'TOTAL'
1334 
1335  ntotp = lx2*ly2*lz2*nelv
1336 
1337  call opscale !normalize vectors to unit norm
1338  $ (vxp(1,jpp),vyp(1,jpp),vzp(1,jpp),tp(1,1,jpp),pertinvnorm)
1339  call cmult(prp(1,jpp),pertinvnorm,ntotp)
1340 
1341  call opscale(exx1p(1,jpp),exy1p(1,jpp),exz1p(1,jpp)
1342  $ ,vgradt1p(1,1,jpp),pertinvnorm)
1343  call opscale(exx2p(1,jpp),exy2p(1,jpp),exz2p(1,jpp)
1344  $ ,vgradt2p(1,1,jpp),pertinvnorm)
1345 
1346  ltotv = lx1*ly1*lz1*lelv
1347  ltotp = lx2*ly2*lz2*lelv
1348 
1349  call cmult( tlagp(1,1,1,jpp),pertinvnorm,ltotv*(lorder-1)*ldimt)
1350  call cmult(vxlagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
1351  call cmult(vylagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
1352  call cmult(vzlagp(1,1,jpp),pertinvnorm,ltotv*(lorder-1))
1353  call cmult(prlagp(1,1,jpp),pertinvnorm,ltotp*(lorder-2))
1354 
1355  if (nio.eq.0) write(6,1) istep,pertnorm,pertinvnorm,jpp,'PNORM'
1356  1 format(i4,1p2e12.4,i4,a5)
1357  pertnorm = pertnorm*pertinvnorm
1358 
1359  return
1360  end
1361 c-----------------------------------------------------------------------
bcdirsc
subroutine bcdirsc(S)
Definition: bdry.f:710
bcneusc
subroutine bcneusc(S, ITYPE)
Definition: bdry.f:786
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
exitt
void exitt()
Definition: comm_mpi.f:604
setqvol
subroutine setqvol(bql)
Definition: conduct.f:131
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
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
make_hpf
subroutine make_hpf
Definition: hpf.f:2
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
opnorm
subroutine opnorm(unorm, ux, uy, uz, type3)
Definition: induct.f:419
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
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
subcol3
subroutine subcol3(a, b, c, n)
Definition: math.f:186
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
add2s1
subroutine add2s1(a, b, c1, n)
Definition: math.f:678
rzero
subroutine rzero(a, n)
Definition: math.f:208
chsign
subroutine chsign(a, n)
Definition: math.f:305
esolver
subroutine esolver(RES, H1, H2, H2INV, INTYPE)
Definition: navier0.f:2
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
opcolv3c
subroutine opcolv3c(a1, a2, a3, b1, b2, b3, c, d)
Definition: navier1.f:2788
convop
subroutine convop(conv, fi)
Definition: navier1.f:3139
opcmult
subroutine opcmult(a, b, c, const)
Definition: navier1.f:2663
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
opcolv3
subroutine opcolv3(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2380
opcolv2
subroutine opcolv2(a1, a2, a3, b1, b2)
Definition: navier1.f:2712
normvc
subroutine normvc(h1, semi, l2, linf, x1, x2, x3)
Definition: navier1.f:2083
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
nekuf
subroutine nekuf(f1, f2, f3)
Definition: navier1.f:1483
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opsub2
subroutine opsub2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2360
opnorm2
function opnorm2(v1, v2, v3)
Definition: pertsupport.f:36
pert_ortho_norm1
subroutine pert_ortho_norm1(k)
Definition: pertsupport.f:550
tnorm
function tnorm(temp)
Definition: pertsupport.f:58
computelyap1
subroutine computelyap1(vxq, vyq, vzq, tq, jpp)
Definition: pertsupport.f:151
dmnorm
real function dmnorm(v1, v2, v3, temp)
Definition: pertsupport.f:70
opscalev
subroutine opscalev(v1, v2, v3, alpha)
Definition: pertsupport.f:120
computelyap
subroutine computelyap
Definition: pertsupport.f:138
rescalepert
subroutine rescalepert(pertnorm, pertinvnorm, jpp)
Definition: pertsupport.f:236
out_pert
subroutine out_pert
Definition: pertsupport.f:454
pert_inner_prod
function pert_inner_prod(i, j)
Definition: pertsupport.f:517
pert_add2s2
subroutine pert_add2s2(i, j, scale)
Definition: pertsupport.f:473
opscale
subroutine opscale(v1, v2, v3, temp, alpha)
Definition: pertsupport.f:100
pert_ortho_norm
subroutine pert_ortho_norm
Definition: pertsupport.f:539
incomprp
subroutine incomprp(ux, uy, uz, up)
Definition: perturb.f:846
cresvipp
subroutine cresvipp(resv1, resv2, resv3, h1, h2)
Definition: perturb.f:511
extrapprp
subroutine extrapprp(prextr)
Definition: perturb.f:927
makeufp
subroutine makeufp
Definition: perturb.f:115
makeabqp
subroutine makeabqp
Definition: perturb.f:751
perturbv
subroutine perturbv(igeom)
Definition: perturb.f:26
lagpresp
subroutine lagpresp
Definition: perturb.f:953
advabp
subroutine advabp
Definition: perturb.f:133
makeqp
subroutine makeqp
Definition: perturb.f:654
makebdqp
subroutine makebdqp
Definition: perturb.f:780
makefp
subroutine makefp
Definition: perturb.f:94
lagfieldp
subroutine lagfieldp
Definition: perturb.f:74
dmnorm1
function dmnorm1(v1, v2, v3, temp)
Definition: perturb.f:1167
convop_adj
subroutine convop_adj(bdux, bduy, bduz, udx, udy, udz, cx, cy, cz)
Definition: perturb.f:318
lagscalp
subroutine lagscalp
Definition: perturb.f:823
heatp
subroutine heatp(igeom)
Definition: perturb.f:544
makextp
subroutine makextp
Definition: perturb.f:428
makeuqp
subroutine makeuqp
Definition: perturb.f:680
fluidp
subroutine fluidp(igeom)
Definition: perturb.f:3
makebdfp
subroutine makebdfp
Definition: perturb.f:467
lyap_scale
subroutine lyap_scale
Definition: perturb.f:968
cdscalp
subroutine cdscalp(igeom)
Definition: perturb.f:579
convabp
subroutine convabp
Definition: perturb.f:705
advabp_adjoint
subroutine advabp_adjoint
Definition: perturb.f:208
outpost2
subroutine outpost2(v1, v2, v3, vp, vt, nfldt, name3)
Definition: prepost.f:1394
settolt
subroutine settolt
Definition: ssolv.f:693
unorm
subroutine unorm
Definition: subs1.f:352
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341