KTH framework for Nek5000 toolboxes; testing version  0.0.1
plan4.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine plan4 (igeom)
3 
4 C Splitting scheme A.G. Tomboulides et al.
5 c Journal of Sci.Comp.,Vol. 12, No. 2, 1998
6 c
7 C NOTE: QTL denotes the so called thermal
8 c divergence and has to be provided
9 c by userqtl.
10 c
11  include 'SIZE'
12  include 'INPUT'
13  include 'GEOM'
14  include 'MASS'
15  include 'SOLN'
16  include 'MVGEOM'
17  include 'TSTEP'
18  include 'ORTHOP'
19  include 'CTIMER'
20 C
21  COMMON /scrns/ res1(lx1,ly1,lz1,lelv)
22  $ , res2(lx1,ly1,lz1,lelv)
23  $ , res3(lx1,ly1,lz1,lelv)
24  $ , dv1(lx1,ly1,lz1,lelv)
25  $ , dv2(lx1,ly1,lz1,lelv)
26  $ , dv3(lx1,ly1,lz1,lelv)
27  $ , respr(lx2,ly2,lz2,lelv)
28  common /scrvh/ h1(lx1,ly1,lz1,lelv)
29  $ , h2(lx1,ly1,lz1,lelv)
30 
31  REAL DPR (LX2,LY2,LZ2,LELV)
32  equivalence(dpr,dv1)
33  LOGICAL IFSTSP
34 
35  REAL DVC (LX1,LY1,LZ1,LELV), DFC(LX1,LY1,LZ1,LELV)
36  REAL DIV1, DIV2, DIF1, DIF2, QTL1, QTL2
37 c
38  intype = -1
39  ntot1 = lx1*ly1*lz1*nelv
40 
41  if (igeom.eq.1) then
42 
43  ! compute explicit contributions bfx,bfy,bfz
44  call makef
45 
46  call sumab(vx_e,vx,vxlag,ntot1,ab,nab)
47  call sumab(vy_e,vy,vylag,ntot1,ab,nab)
48  if (if3d) call sumab(vz_e,vz,vzlag,ntot1,ab,nab)
49 
50  else
51 
52  if(iflomach) call opcolv(bfx,bfy,bfz,vtrans)
53 
54  ! add user defined divergence to qtl
55  call add2 (qtl,usrdiv,ntot1)
56 
57  if (igeom.eq.2) call lagvel
58 
59  ! mask Dirichlet boundaries
60  call bcdirvc (vx,vy,vz,v1mask,v2mask,v3mask)
61 
62  ! compute pressure
63  call copy(prlag,pr,ntot1)
64  if (icalld.eq.0) tpres=0.0
65  icalld=icalld+1
66  npres=icalld
67  etime1=dnekclock()
68 
69  call crespsp (respr)
70  call invers2 (h1,vtrans,ntot1)
71  call rzero (h2,ntot1)
72  call ctolspl (tolspl,respr)
73  napproxp(1) = laxtp
74  call hsolve ('PRES',dpr,respr,h1,h2
75  $ ,pmask,vmult
76  $ ,imesh,tolspl,nmxp,1
77  $ ,approxp,napproxp,binvm1)
78  call add2 (pr,dpr,ntot1)
79  call ortho (pr)
80 
81  tpres=tpres+(dnekclock()-etime1)
82 
83  ! compute velocity
84  if(ifstrs .and. .not.ifaxis) then
85  call bcneutr
86  call cresvsp_weak(res1,res2,res3,h1,h2)
87  else
88  call cresvsp (res1,res2,res3,h1,h2)
89  endif
90  call ophinv (dv1,dv2,dv3,res1,res2,res3,h1,h2,tolhv,nmxv)
91  call opadd2 (vx,vy,vz,dv1,dv2,dv3)
92 
93  endif
94 
95  return
96  END
97 
98 c-----------------------------------------------------------------------
99  subroutine crespsp (respr)
100 
101 C Compute startresidual/right-hand-side in the pressure
102 
103  include 'SIZE'
104  include 'TOTAL'
105 
106  REAL RESPR (LX1*LY1*LZ1,LELV)
107 c
108  COMMON /scrns/ ta1(lx1*ly1*lz1,lelv)
109  $ , ta2(lx1*ly1*lz1,lelv)
110  $ , ta3(lx1*ly1*lz1,lelv)
111  $ , wa1(lx1*ly1*lz1*lelv)
112  $ , wa2(lx1*ly1*lz1*lelv)
113  $ , wa3(lx1*ly1*lz1*lelv)
114  COMMON /scrmg/ w1(lx1*ly1*lz1,lelv)
115  $ , w2(lx1*ly1*lz1,lelv)
116  $ , w3(lx1*ly1*lz1,lelv)
117 
118  common /scruz/ sij(lx1*ly1*lz1,6,lelv)
119  parameter(lr=lx1*ly1*lz1)
120  common /scrvz/ ur(lr),us(lr),ut(lr)
121  $ , vr(lr),vs(lr),vt(lr)
122  $ , wr(lr),ws(lr),wt(lr)
123 
124  CHARACTER CB*3
125 
126  nxyz1 = lx1*ly1*lz1
127  ntot1 = nxyz1*nelv
128  nfaces = 2*ldim
129 
130 c -mu*curl(curl(v))
131  call op_curl (ta1,ta2,ta3,vx_e,vy_e,vz_e,
132  & .true.,w1,w2)
133  if(ifaxis) then
134  CALL col2 (ta2, omask,ntot1)
135  CALL col2 (ta3, omask,ntot1)
136  endif
137  call op_curl (wa1,wa2,wa3,ta1,ta2,ta3,.true.,w1,w2)
138  if(ifaxis) then
139  CALL col2 (wa2, omask,ntot1)
140  CALL col2 (wa3, omask,ntot1)
141  endif
142  call opcolv (wa1,wa2,wa3,bm1)
143 c
144  call opgrad (ta1,ta2,ta3,qtl)
145  if(ifaxis) then
146  CALL col2 (ta2, omask,ntot1)
147  CALL col2 (ta3, omask,ntot1)
148  endif
149  scale = -4./3.
150  call opadd2cm (wa1,wa2,wa3,ta1,ta2,ta3,scale)
151 
152 c compute stress tensor for ifstrs formulation - variable viscosity Pn-Pn
153  if (ifstrs .and. ifvvisp) then
154  call opgrad (ta1,ta2,ta3,vdiff)
155  call invcol2 (ta1,vdiff,ntot1)
156  call invcol2 (ta2,vdiff,ntot1)
157  call invcol2 (ta3,vdiff,ntot1)
158 
159  nij = 3
160  if (if3d.or.ifaxis) nij=6
161 
162  call comp_sij (sij,nij,vx_e,vy_e,vz_e,
163  & ur,us,ut,vr,vs,vt,wr,ws,wt)
164  call col_mu_sij (w1,w2,w3,ta1,ta2,ta3,sij,nij)
165 
166  call opcolv (ta1,ta2,ta3,qtl)
167  scale2 = -2./3.
168  call opadd2cm (w1,w2,w3,ta1,ta2,ta3,scale2)
169  call opsub2 (wa1,wa2,wa3,w1,w2,w3)
170 
171  endif
172 
173  call invcol3 (w1,vdiff,vtrans,ntot1)
174  call opcolv (wa1,wa2,wa3,w1)
175 
176 c add old pressure term because we solve for delta p
177  call invers2 (ta1,vtrans,ntot1)
178  call rzero (ta2,ntot1)
179 
180  call bcdirpr (pr)
181 
182  call axhelm (respr,pr,ta1,ta2,imesh,1)
183  call chsign (respr,ntot1)
184 
185 c add explicit (NONLINEAR) terms
186  n = lx1*ly1*lz1*nelv
187  do i=1,n
188  ta1(i,1) = bfx(i,1,1,1)/vtrans(i,1,1,1,1)-wa1(i)
189  ta2(i,1) = bfy(i,1,1,1)/vtrans(i,1,1,1,1)-wa2(i)
190  ta3(i,1) = bfz(i,1,1,1)/vtrans(i,1,1,1,1)-wa3(i)
191  enddo
192  call opdssum (ta1,ta2,ta3)
193  do i=1,n
194  ta1(i,1) = ta1(i,1)*binvm1(i,1,1,1)
195  ta2(i,1) = ta2(i,1)*binvm1(i,1,1,1)
196  ta3(i,1) = ta3(i,1)*binvm1(i,1,1,1)
197  enddo
198 
199  if (if3d) then
200  call cdtp (wa1,ta1,rxm1,sxm1,txm1,1)
201  call cdtp (wa2,ta2,rym1,sym1,tym1,1)
202  call cdtp (wa3,ta3,rzm1,szm1,tzm1,1)
203  do i=1,n
204  respr(i,1) = respr(i,1)+wa1(i)+wa2(i)+wa3(i)
205  enddo
206  else
207  call cdtp (wa1,ta1,rxm1,sxm1,txm1,1)
208  call cdtp (wa2,ta2,rym1,sym1,tym1,1)
209 
210  do i=1,n
211  respr(i,1) = respr(i,1)+wa1(i)+wa2(i)
212  enddo
213  endif
214 
215 C add thermal divergence
216  dtbd = bd(1)/dt
217  call admcol3(respr,qtl,bm1,dtbd,ntot1)
218 
219 C surface terms
220  DO 100 iel=1,nelv
221  DO 300 ifc=1,nfaces
222  CALL rzero (w1(1,iel),nxyz1)
223  CALL rzero (w2(1,iel),nxyz1)
224  IF (ldim.EQ.3)
225  $ CALL rzero (w3(1,iel),nxyz1)
226  cb = cbc(ifc,iel,ifield)
227  IF (cb(1:1).EQ.'V'.OR.cb(1:1).EQ.'v'.or.
228  $ cb.eq.'MV '.or.cb.eq.'mv '.or.cb.eq.'shl') then
229  CALL faccl3
230  $ (w1(1,iel),vx(1,1,1,iel),unx(1,1,ifc,iel),ifc)
231  CALL faccl3
232  $ (w2(1,iel),vy(1,1,1,iel),uny(1,1,ifc,iel),ifc)
233  IF (ldim.EQ.3)
234  $ CALL faccl3
235  $ (w3(1,iel),vz(1,1,1,iel),unz(1,1,ifc,iel),ifc)
236  ELSE IF (cb(1:3).EQ.'SYM') THEN
237  CALL faccl3
238  $ (w1(1,iel),ta1(1,iel),unx(1,1,ifc,iel),ifc)
239  CALL faccl3
240  $ (w2(1,iel),ta2(1,iel),uny(1,1,ifc,iel),ifc)
241  IF (ldim.EQ.3)
242  $ CALL faccl3
243  $ (w3(1,iel),ta3(1,iel),unz(1,1,ifc,iel),ifc)
244  ENDIF
245  CALL add2 (w1(1,iel),w2(1,iel),nxyz1)
246  IF (ldim.EQ.3)
247  $ CALL add2 (w1(1,iel),w3(1,iel),nxyz1)
248  CALL faccl2 (w1(1,iel),area(1,1,ifc,iel),ifc)
249  IF (cb(1:1).EQ.'V'.OR.cb(1:1).EQ.'v'.or.
250  $ cb.eq.'MV '.or.cb.eq.'mv ') then
251  CALL cmult(w1(1,iel),dtbd,nxyz1)
252  endif
253  CALL sub2 (respr(1,iel),w1(1,iel),nxyz1)
254  300 CONTINUE
255  100 CONTINUE
256 
257 C Assure that the residual is orthogonal to (1,1,...,1)T
258 C (only if all Dirichlet b.c.)
259  CALL ortho (respr)
260 
261  return
262  END
263 c----------------------------------------------------------------------
264  subroutine cresvsp (resv1,resv2,resv3,h1,h2)
265 
266 C Compute the residual for the velocity
267 
268  include 'SIZE'
269  include 'TOTAL'
270 
271  real resv1(lx1,ly1,lz1,lelv)
272  $ , resv2(lx1,ly1,lz1,lelv)
273  $ , resv3(lx1,ly1,lz1,lelv)
274  $ , h1 (lx1,ly1,lz1,lelv)
275  $ , h2 (lx1,ly1,lz1,lelv)
276 
277  COMMON /scruz/ ta1(lx1,ly1,lz1,lelv)
278  $ , ta2(lx1,ly1,lz1,lelv)
279  $ , ta3(lx1,ly1,lz1,lelv)
280  $ , ta4(lx1,ly1,lz1,lelv)
281 
282  ntot = lx1*ly1*lz1*nelv
283  intype = -1
284 
285  CALL sethlm (h1,h2,intype)
286 
287  CALL ophx (resv1,resv2,resv3,vx,vy,vz,h1,h2)
288  CALL opchsgn (resv1,resv2,resv3)
289 
290  scale = -1./3.
291  if (ifstrs) scale = 2./3.
292 
293  call col3 (ta4,vdiff,qtl,ntot)
294  call add2s1 (ta4,pr,scale,ntot)
295  call opgrad (ta1,ta2,ta3,ta4)
296  if(ifaxis) then
297  CALL col2 (ta2, omask,ntot)
298  CALL col2 (ta3, omask,ntot)
299  endif
300 c
301  call opsub2 (resv1,resv2,resv3,ta1,ta2,ta3)
302  call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
303 
304  return
305  end
306 
307 c----------------------------------------------------------------------
308  subroutine cresvsp_weak (resv1,resv2,resv3,h1,h2)
309 
310 C Compute the residual for the velocity
311 
312  include 'SIZE'
313  include 'TOTAL'
314 
315  real resv1(lx1,ly1,lz1,lelv)
316  $ , resv2(lx1,ly1,lz1,lelv)
317  $ , resv3(lx1,ly1,lz1,lelv)
318  $ , h1 (lx1,ly1,lz1,lelv)
319  $ , h2 (lx1,ly1,lz1,lelv)
320 
321  COMMON /scruz/ ta1(lx1,ly1,lz1,lelv)
322  $ , ta2(lx1,ly1,lz1,lelv)
323  $ , ta3(lx1,ly1,lz1,lelv)
324  $ , ta4(lx1,ly1,lz1,lelv)
325  COMMON /scrmg/ wa1(lx1*ly1*lz1,lelv)
326  $ , wa2(lx1*ly1*lz1,lelv)
327  $ , wa3(lx1*ly1*lz1,lelv)
328 
329  ntot = lx1*ly1*lz1*nelv
330  intype = -1
331 
332  CALL oprzero (resv1,resv2,resv3)
333  CALL oprzero (wa1 ,wa2 ,wa3 )
334  CALL oprzero (ta1 ,ta2 ,ta3 )
335 
336  CALL sethlm (h1,h2,intype)
337 
338  CALL ophx (resv1,resv2,resv3,vx,vy,vz,h1,h2)
339  CALL opchsgn (resv1,resv2,resv3)
340 
341  scale = -1./3.
342  if (ifstrs) scale = 2./3.
343 
344  call col3 (ta4,vdiff,qtl,ntot)
345  call cmult (ta4,scale,ntot)
346  call opgrad (ta1,ta2,ta3,ta4)
347 
348  call cdtp (wa1,pr ,rxm1,sxm1,txm1,1)
349  call cdtp (wa2,pr ,rym1,sym1,tym1,1)
350  if(if3d) call cdtp (wa3,pr ,rzm1,szm1,tzm1,1)
351 
352  call sub2 (ta1,wa1,ntot)
353  call sub2 (ta2,wa2,ntot)
354  if(if3d) call sub2 (ta3,wa3,ntot)
355 
356  if(ifaxis) then
357  CALL col2 (ta2, omask,ntot)
358  CALL col2 (ta3, omask,ntot)
359  endif
360 c
361  call opsub2 (resv1,resv2,resv3,ta1,ta2,ta3)
362  call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
363 
364  return
365  end
366 
367 c-----------------------------------------------------------------------
368  subroutine op_curl(w1,w2,w3,u1,u2,u3,ifavg,work1,work2)
369 c
370  include 'SIZE'
371  include 'TOTAL'
372 c
373  real duax(lx1), ta(lx1,ly1,lz1,lelv)
374 
375  logical ifavg
376 c
377  real w1(1),w2(1),w3(1),work1(1),work2(1),u1(1),u2(1),u3(1)
378 c
379  ntot = lx1*ly1*lz1*nelv
380  nxyz = lx1*ly1*lz1
381 c work1=dw/dy ; work2=dv/dz
382  call dudxyz(work1,u3,rym1,sym1,tym1,jacm1,1,2)
383  if (if3d) then
384  call dudxyz(work2,u2,rzm1,szm1,tzm1,jacm1,1,3)
385  call sub3(w1,work1,work2,ntot)
386  else
387  call copy(w1,work1,ntot)
388 
389  if(ifaxis) then
390  call copy (ta,u3,ntot)
391  do iel = 1,nelv
392  if(ifrzer(iel)) then
393  call rzero (ta(1,1,1,iel),lx1)
394  call mxm (ta(1,1,1,iel),lx1,datm1,ly1,duax,1)
395  call copy (ta(1,1,1,iel),duax,lx1)
396  endif
397  call col2 (ta(1,1,1,iel),yinvm1(1,1,1,iel),nxyz)
398  enddo
399  call add2 (w1,ta,ntot)
400  endif
401 
402  endif
403 c work1=du/dz ; work2=dw/dx
404  if (if3d) then
405  call dudxyz(work1,u1,rzm1,szm1,tzm1,jacm1,1,3)
406  call dudxyz(work2,u3,rxm1,sxm1,txm1,jacm1,1,1)
407  call sub3(w2,work1,work2,ntot)
408  else
409  call rzero (work1,ntot)
410  call dudxyz(work2,u3,rxm1,sxm1,txm1,jacm1,1,1)
411  call sub3(w2,work1,work2,ntot)
412  endif
413 c work1=dv/dx ; work2=du/dy
414  call dudxyz(work1,u2,rxm1,sxm1,txm1,jacm1,1,1)
415  call dudxyz(work2,u1,rym1,sym1,tym1,jacm1,1,2)
416  call sub3(w3,work1,work2,ntot)
417 c
418 c Avg at bndry
419 c
420 c if (ifavg) then
421  if (ifavg .and. .not. ifcyclic) then
422 
423  ifielt = ifield
424  ifield = 1
425 
426  call opcolv (w1,w2,w3,bm1)
427  call opdssum (w1,w2,w3)
428  call opcolv (w1,w2,w3,binvm1)
429 
430  ifield = ifielt
431 
432  endif
433 c
434  return
435  end
436 
437 c-----------------------------------------------------------------------
438  subroutine opadd2cm (a1,a2,a3,b1,b2,b3,c)
439  include 'SIZE'
440  REAL A1(1),A2(1),A3(1),B1(1),B2(1),B3(1),C
441  ntot1=lx1*ly1*lz1*nelv
442  if (ldim.eq.3) then
443  do i=1,ntot1
444  a1(i) = a1(i) + b1(i)*c
445  a2(i) = a2(i) + b2(i)*c
446  a3(i) = a3(i) + b3(i)*c
447  enddo
448  else
449  do i=1,ntot1
450  a1(i) = a1(i) + b1(i)*c
451  a2(i) = a2(i) + b2(i)*c
452  enddo
453  endif
454  return
455  END
456 
457 c-----------------------------------------------------------------------
458  subroutine split_vis
459 
460 c Split viscosity into a constant implicit (VDIFF) and variable
461 c explicit (VDIFF_E) part.
462 
463  include 'SIZE'
464  include 'TOTAL'
465 
466  n = lx1*ly1*lz1*nelv
467 
468  dnu_star = -nu_star
469  call cadd2 (vdiff_e,vdiff,dnu_star,n) ! set explicit part
470 
471  call cfill (vdiff,nu_star,n) ! set implicit part
472 
473  return
474  end
475 c-----------------------------------------------------------------------
476  subroutine redo_split_vis ! Redo split viscosity
477 
478  include 'SIZE'
479  include 'TOTAL'
480 
481  n = lx1*ly1*lz1*nelv
482  call add2(vdiff,vdiff_e,n) ! sum up explicit and implicit part
483 
484  return
485  end
486 c-----------------------------------------------------------------------
487  subroutine col_mu_sij(w1,w2,w3,ta1,ta2,ta3,sij,nij)
488 
489  include 'SIZE'
490  include 'TOTAL'
491 
492  parameter(lr=lx1*ly1*lz1)
493  real w1 (lr,1),w2 (lr,1),w3 (lr,1)
494  real ta1(lr,1),ta2(lr,1),ta3(lr,1)
495  real sij(lr,nij,1)
496  integer e
497 
498  nxyz1 = lx1*ly1*lz1
499 
500  if (if3d.or.ifaxis) then
501  do e=1,nelv
502  call vdot3 (w1(1,e),ta1(1,e),ta2(1,e),ta3(1,e)
503  $ ,sij(1,1,e),sij(1,4,e),sij(1,6,e),nxyz1)
504  call vdot3 (w2(1,e),ta1(1,e),ta2(1,e),ta3(1,e)
505  $ ,sij(1,4,e),sij(1,2,e),sij(1,5,e),nxyz1)
506  call vdot3 (w3(1,e),ta1(1,e),ta2(1,e),ta3(1,e)
507  $ ,sij(1,6,e),sij(1,5,e),sij(1,3,e),nxyz1)
508  enddo
509 
510  else
511 
512  do e=1,nelv
513  call vdot2 (w1(1,e),ta1(1,e),ta2(1,e)
514  $ ,sij(1,1,e),sij(1,3,e),nxyz1)
515  call vdot2 (w2,ta1(1,e),ta2(1,e)
516  $ ,sij(1,3,e),sij(1,2,e),nxyz1)
517  call rzero (w3(1,e),nxyz1)
518  enddo
519 
520  endif
521 
522  return
523  end
524 c-----------------------------------------------------------------------
525  subroutine sumab(v,vv,vvlag,ntot,ab_,nab_)
526 c
527 c sum up AB/BDF contributions
528 c
529  include 'SIZE'
530 
531  real vvlag(lx1*ly1*lz1*lelv,*)
532  real ab_(*)
533 
534  ab0 = ab_(1)
535  ab1 = ab_(2)
536  ab2 = ab_(3)
537 
538  call add3s2(v,vv,vvlag(1,1),ab0,ab1,ntot)
539  if(nab_.eq.3) call add2s2(v,vvlag(1,2),ab2,ntot)
540 
541  return
542  end
543 c-----------------------------------------------------------------------
544  subroutine userqtl_scig
545 c
546 c Compute the thermal divergence QTL for an ideal single component gas
547 c QTL := div(v) = -1/rho * Drho/Dt
548 c
549  include 'SIZE'
550  include 'TOTAL'
551  include 'CVODE'
552 
553  common /scrns/ w1(lx1,ly1,lz1,lelt)
554  $ ,w2(lx1,ly1,lz1,lelt)
555  $ ,w3(lx1,ly1,lz1,lelt)
556  $ ,tx(lx1,ly1,lz1,lelt)
557  $ ,ty(lx1,ly1,lz1,lelt)
558  $ ,tz(lx1,ly1,lz1,lelt)
559 
560  nxyz = lx1*ly1*lz1
561  ntot = nxyz*nelv
562 
563  ifld_save = ifield
564 
565 c - - Assemble RHS of T-eqn
566  ifield=2
567 
568  call makeuq
569  call copy(qtl,bq,ntot)
570 
571  ifield=1 !set right gs handle (QTL is only defined on the velocity mesh)
572  call opgrad (tx,ty,tz,t)
573  call opdssum (tx,ty,tz)
574  call opcolv (tx,ty,tz,binvm1)
575  call opcolv (tx,ty,tz,vdiff(1,1,1,1,2))
576  call opdiv (w2,tx,ty,tz)
577 
578  call add2 (qtl,w2,ntot)
579  call dssum (qtl,lx1,ly1,lz1)
580  call col2 (qtl,binvm1,ntot)
581 
582  ! QTL = T_RHS/(rho*cp**T)
583  call col3 (w2,vtrans(1,1,1,1,2),t,ntot)
584  call invcol2 (qtl,w2,ntot)
585 
586  dp0thdt = 0.0
587  if (ifdp0dt) then
588  dd = (1.0 - gamma0)/gamma0
589  call rone(w1,ntot)
590  call cmult(w1,dd,ntot)
591 
592  call invcol3(w2,vtrans(1,1,1,1,2),vtrans,ntot)
593  call invcol2(w1,w2,ntot)
594 
595  call cadd(w1,1.0,ntot)
596  call copy(w2,w1,ntot)
597  call col2(w1,bm1,ntot)
598 
599  p0alph1 = 1. / glsum(w1,ntot)
600 
601  call copy (w1,qtl,ntot)
602  call col2 (w1,bm1,ntot)
603 
604  termq = glsum(w1,ntot)
605  if (ifcvfun) then
606  termv = glcflux(vx,vy,vz)
607  prhs = p0alph1*(termq - termv)
608  pcoef =(cv_bd(1) - cv_dtnek*prhs)
609  call add3s2(saqpq,p0thn,p0thlag(1),cv_bd(2),cv_bd(3),1)
610  if(nbd.eq.3) call add2s2(saqpq,p0thlag(2),cv_bd(4),1)
611  p0th = saqpq / pcoef
612  else
613  termv = glcflux(vx_e,vy_e,vz_e)
614  prhs = p0alph1*(termq - termv)
615  pcoef =(bd(1) - dt*prhs)
616  call add3s2(saqpq,p0thn,p0thlag(1),bd(2),bd(3),1)
617  if(nbd.eq.3) call add2s2(saqpq,p0thlag(2),bd(4),1)
618  p0th = saqpq / pcoef
619  p0thlag(2) = p0thlag(1)
620  p0thlag(1) = p0thn
621  p0thn = p0th
622  endif
623 
624  dp0thdt= prhs*p0th
625 
626  dd =-prhs
627  call cmult(w2,dd,ntot)
628  call add2 (qtl,w2,ntot)
629  endif
630 
631  ifield = ifld_save
632 
633  return
634  end
635 c-----------------------------------------------------------------------
636  subroutine qthermal
637 c
638 c generic qtl wrapper
639 c
640  include 'SIZE'
641  include 'TOTAL'
642 
643  ntot = lx1*ly1*lz1*nelv
644 
645  call rzero(qtl,ntot)
646  call userqtl()
647 
648  return
649  end
650 c-----------------------------------------------------------------------
651  subroutine printdiverr
652 c
653  include 'SIZE'
654  include 'TOTAL'
655 
656  COMMON /scrns/ dvc(lx1,ly1,lz1,lelv),
657  $ dv1(lx1,ly1,lz1,lelv),
658  $ dv2(lx1,ly1,lz1,lelv),
659  $ dfc(lx1,ly1,lz1,lelv)
660 
661  ntot1 = lx1*ly1*lz1*nelv
662 
663  !Calculate Divergence norms of new VX,VY,VZ
664  CALL opdiv (dvc,vx,vy,vz)
665  CALL dssum (dvc,lx1,ly1,lz1)
666  CALL col2 (dvc,binvm1,ntot1)
667 
668  CALL col3 (dv1,dvc,bm1,ntot1)
669  div1 = glsum(dv1,ntot1)/volvm1
670 
671  CALL col3 (dv2,dvc,dvc,ntot1)
672  CALL col2 (dv2,bm1 ,ntot1)
673  div2 = glsum(dv2,ntot1)/volvm1
674  div2 = sqrt(div2)
675 
676  !Calculate Divergence difference norms
677  CALL sub3 (dfc,dvc,qtl,ntot1)
678  CALL col3 (dv1,dfc,bm1,ntot1)
679  dif1 = glsum(dv1,ntot1)/volvm1
680 
681  CALL col3 (dv2,dfc,dfc,ntot1)
682  CALL col2 (dv2,bm1 ,ntot1)
683  dif2 = glsum(dv2,ntot1)/volvm1
684  dif2 = sqrt(dif2)
685 
686  CALL col3 (dv1,qtl,bm1,ntot1)
687  qtl1 = glsum(dv1,ntot1)/volvm1
688 
689  CALL col3 (dv2,qtl,qtl,ntot1)
690  CALL col2 (dv2,bm1 ,ntot1)
691  qtl2 = glsum(dv2,ntot1)/volvm1
692  qtl2 = sqrt(qtl2)
693 
694  IF (nio.EQ.0) THEN
695  WRITE(6,'(13X,A,1p2e13.4)')
696  & 'L1/L2 DIV(V) ',div1,div2
697  WRITE(6,'(13X,A,1p2e13.4)')
698  & 'L1/L2 QTL ',qtl1,qtl2
699  WRITE(6,'(13X,A,1p2e13.4)')
700  & 'L1/L2 DIV(V)-QTL ',dif1,dif2
701  IF (dif2.GT.0.1) WRITE(6,'(13X,A)')
702  & 'WARNING: DIV(V)-QTL too large!'
703  ENDIF
704 
705  return
706  end
707 c-----------------------------------------------------------------------
708  SUBROUTINE bcdirpr(S)
709 C
710 C Apply Dirichlet boundary conditions to surface of Pressure.
711 C Use IFIELD=1.
712 C
713  include 'SIZE'
714  include 'TSTEP'
715  include 'INPUT'
716  include 'SOLN'
717  include 'TOPOL'
718  include 'CTIMER'
719 C
720  dimension s(lx1,ly1,lz1,lelt)
721  COMMON /scrsf/ tmp(lx1,ly1,lz1,lelt)
722  $ , tma(lx1,ly1,lz1,lelt)
723  $ , smu(lx1,ly1,lz1,lelt)
724  common /nekcb/ cb
725  CHARACTER CB*3
726 
727  if (icalld.eq.0) then
728  tusbc=0.0
729  nusbc=0
730  icalld=icalld+1
731  endif
732  nusbc=nusbc+1
733  etime1=dnekclock()
734 C
735  ifld = 1
736  nfaces = 2*ldim
737  nxyz = lx1*ly1*lz1
738  nel = nelfld(ifield)
739  ntot = nxyz*nel
740  nfldt = nfield - 1
741 C
742  CALL rzero(tmp,ntot)
743 C
744 C pressure boundary condition
745 C
746  DO 2100 isweep=1,2
747 C
748  DO 2010 ie=1,nel
749  DO 2010 iface=1,nfaces
750  cb=cbc(iface,ie,ifield)
751  bc1=bc(1,iface,ie,ifield)
752  IF (cb.EQ.'O ' .or. cb.eq.'ON ' .or.
753  $ cb.eq.'o ' .or. cb.eq.'on ')
754  $ CALL faceis (cb,tmp(1,1,1,ie),ie,iface,lx1,ly1,lz1)
755  2010 CONTINUE
756 C
757 C Take care of Neumann-Dirichlet shared edges...
758 C
759  IF (isweep.EQ.1) CALL dsop(tmp,'MXA',lx1,ly1,lz1)
760  IF (isweep.EQ.2) CALL dsop(tmp,'MNA',lx1,ly1,lz1)
761  2100 CONTINUE
762 C
763 C Copy temporary array to temperature array.
764 C
765  CALL col2(s,pmask,ntot)
766  CALL add2(s,tmp,ntot)
767 
768  tusbc=tusbc+(dnekclock()-etime1)
769 
770  RETURN
771  END
772 C
773 c-----------------------------------------------------------------------
faceis
subroutine faceis(CB, S, IEL, IFACE, NX, NY, NZ)
Definition: bdry.f:919
bcdirvc
subroutine bcdirvc(V1, V2, V3, mask1, mask2, mask3)
Definition: bdry.f:574
glcflux
real function glcflux(tx, ty, tz)
Definition: bdry.f:2015
bcneutr
subroutine bcneutr
Definition: bdry.f:1200
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
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
cadd
subroutine cadd(a, const, n)
Definition: math.f:326
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
admcol3
subroutine admcol3(a, b, c, d, n)
Definition: math.f:1556
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
vdot2
subroutine vdot2(dot, u1, u2, v1, v2, n)
Definition: math.f:447
cadd2
subroutine cadd2(a, b, const, n)
Definition: math.f:346
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
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
add2s1
subroutine add2s1(a, b, c1, n)
Definition: math.f:678
vdot3
subroutine vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Definition: math.f:462
rzero
subroutine rzero(a, n)
Definition: math.f:208
chsign
subroutine chsign(a, n)
Definition: math.f:305
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
cdtp
subroutine cdtp(dtx, x, rm2, sm2, tm2, isd)
Definition: navier1.f:331
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
ortho
subroutine ortho(respr)
Definition: navier1.f:224
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
makef
subroutine makef
Definition: navier1.f:1426
lagvel
subroutine lagvel
Definition: navier1.f:1930
opcolv
subroutine opcolv(a1, a2, a3, c)
Definition: navier1.f:2418
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
ctolspl
subroutine ctolspl(tolspl, respr)
Definition: navier1.f:194
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
opchsgn
subroutine opchsgn(a, b, c)
Definition: navier1.f:2464
opsub2
subroutine opsub2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2360
oprzero
subroutine oprzero(a, b, c)
Definition: navier1.f:2643
opgrad
subroutine opgrad(out1, out2, out3, inp)
Definition: navier1.f:298
hsolve
subroutine hsolve(name, u, r, h1, h2, vmk, vml, imsh, tol, maxit, isd, approx, napprox, bi)
Definition: navier4.f:564
comp_sij
subroutine comp_sij(sij, nij, u, v, w, ur, us, ut, vr, vs, vt, wr, ws, wt)
Definition: navier5.f:1824
crespsp
subroutine crespsp(respr)
Definition: plan4.f:100
op_curl
subroutine op_curl(w1, w2, w3, u1, u2, u3, ifavg, work1, work2)
Definition: plan4.f:369
plan4
subroutine plan4(igeom)
Definition: plan4.f:3
qthermal
subroutine qthermal
Definition: plan4.f:637
userqtl_scig
subroutine userqtl_scig
Definition: plan4.f:545
split_vis
subroutine split_vis
Definition: plan4.f:459
cresvsp_weak
subroutine cresvsp_weak(resv1, resv2, resv3, h1, h2)
Definition: plan4.f:309
redo_split_vis
subroutine redo_split_vis
Definition: plan4.f:477
printdiverr
subroutine printdiverr
Definition: plan4.f:652
sumab
subroutine sumab(v, vv, vvlag, ntot, ab_, nab_)
Definition: plan4.f:526
opadd2cm
subroutine opadd2cm(a1, a2, a3, b1, b2, b3, c)
Definition: plan4.f:439
cresvsp
subroutine cresvsp(resv1, resv2, resv3, h1, h2)
Definition: plan4.f:265
col_mu_sij
subroutine col_mu_sij(w1, w2, w3, ta1, ta2, ta3, sij, nij)
Definition: plan4.f:488
bcdirpr
subroutine bcdirpr(S)
Definition: plan4.f:709
faccl2
subroutine faccl2(a, b, iface1)
Definition: subs1.f:909
faccl3
subroutine faccl3(a, b, c, iface1)
Definition: subs1.f:944
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014