source: CIVL/examples/fortran/nek5000/core/plan4.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 21.4 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine plan4 (igeom)
3
4C Splitting scheme A.G. Tomboulides et al.
5c Journal of Sci.Comp.,Vol. 12, No. 2, 1998
6c
7C NOTE: QTL denotes the so called thermal
8c divergence and has to be provided
9c by userqtl.
10c
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'
20C
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
37c
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
98c-----------------------------------------------------------------------
99 subroutine crespsp (respr)
100
101C Compute startresidual/right-hand-side in the pressure
102
103 INCLUDE 'SIZE'
104 INCLUDE 'TOTAL'
105
106 REAL RESPR (LX1*LY1*LZ1,LELV)
107c
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
130c -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)
143c
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
152c 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
176c 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
185c 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
215C add thermal divergence
216 dtbd = BD(1)/DT
217 call admcol3(respr,QTL,bm1,dtbd,ntot1)
218
219C 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
257C Assure that the residual is orthogonal to (1,1,...,1)T
258C (only if all Dirichlet b.c.)
259 CALL ORTHO (RESPR)
260
261 return
262 END
263c----------------------------------------------------------------------
264 subroutine cresvsp (resv1,resv2,resv3,h1,h2)
265
266C 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
300c
301 call opsub2 (resv1,resv2,resv3,ta1,ta2,ta3)
302 call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
303
304 return
305 end
306
307c----------------------------------------------------------------------
308 subroutine cresvsp_weak (resv1,resv2,resv3,h1,h2)
309
310C 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
360c
361 call opsub2 (resv1,resv2,resv3,ta1,ta2,ta3)
362 call opadd2 (resv1,resv2,resv3,bfx,bfy,bfz)
363
364 return
365 end
366
367c-----------------------------------------------------------------------
368 subroutine op_curl(w1,w2,w3,u1,u2,u3,ifavg,work1,work2)
369c
370 include 'SIZE'
371 include 'TOTAL'
372c
373 real duax(lx1), ta(lx1,ly1,lz1,lelv)
374
375 logical ifavg
376c
377 real w1(1),w2(1),w3(1),work1(1),work2(1),u1(1),u2(1),u3(1)
378c
379 ntot = lx1*ly1*lz1*nelv
380 nxyz = lx1*ly1*lz1
381c 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
403c 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
413c 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)
417c
418c Avg at bndry
419c
420c 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
433c
434 return
435 end
436
437c-----------------------------------------------------------------------
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
457c-----------------------------------------------------------------------
458 subroutine split_vis
459
460c Split viscosity into a constant implicit (VDIFF) and variable
461c 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
475c-----------------------------------------------------------------------
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
486c-----------------------------------------------------------------------
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
524c-----------------------------------------------------------------------
525 subroutine sumab(v,vv,vvlag,ntot,ab_,nab_)
526c
527c sum up AB/BDF contributions
528c
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
543c-----------------------------------------------------------------------
544 subroutine userqtl_scig
545c
546c Compute the thermal divergence QTL for an ideal single component gas
547c QTL := div(v) = -1/rho * Drho/Dt
548c
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
565c - - 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
635c-----------------------------------------------------------------------
636 subroutine qthermal
637c
638c generic qtl wrapper
639c
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
650c-----------------------------------------------------------------------
651 subroutine printdiverr
652c
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
707c-----------------------------------------------------------------------
708 SUBROUTINE BCDIRPR(S)
709C
710C Apply Dirichlet boundary conditions to surface of Pressure.
711C Use IFIELD=1.
712C
713 INCLUDE 'SIZE'
714 INCLUDE 'TSTEP'
715 INCLUDE 'INPUT'
716 INCLUDE 'SOLN'
717 INCLUDE 'TOPOL'
718 INCLUDE 'CTIMER'
719C
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()
734C
735 IFLD = 1
736 NFACES = 2*ldim
737 NXYZ = lx1*ly1*lz1
738 NEL = NELFLD(IFIELD)
739 NTOT = NXYZ*NEL
740 NFLDT = NFIELD - 1
741C
742 CALL RZERO(TMP,NTOT)
743C
744C pressure boundary condition
745C
746 DO 2100 ISWEEP=1,2
747C
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
756C
757C Take care of Neumann-Dirichlet shared edges...
758C
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
762C
763C Copy temporary array to temperature array.
764C
765 CALL COL2(S,PMASK,NTOT)
766 CALL ADD2(S,TMP,NTOT)
767
768 tusbc=tusbc+(dnekclock()-etime1)
769
770 RETURN
771 END
772C
773c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.