Actual source code: user.F
1: !---------------------------------------------------------------
2: ! The following subroutines are from node2t.f in the original
3: ! code - D. K. Kaushik (1/17/97)
4: !---------------------------------------------------------------
5: !
6: !
7: ! 2-D Navier Stokes on Unstructured Grids
9: !============================== FORLINK ==============================72
10: !
11: ! FORLINK establishes links between FORTRAN common blocks and C
12: !
13: !=====================================================================72
14: subroutine FORLINK()
16: implicit none
17: #include <petsc/finclude/petscsys.h>
18: PetscScalar title(20),beta,alpha
19: PetscScalar Re,dt,tot,res0,resc
20: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
21: PetscScalar cfl1,cfl2
22: integer nsmoth,iflim,itran,nbtran,jupdate, &
23: & nstage,ncyct,iramp,nitfo
24: PetscScalar gtol
25: integer icycle,nsrch,ilu0,ifcn
26: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
27: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
28: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
29: & nstage,ncyct,iramp,nitfo
30: common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
32: call CLINK(title,cfl1,gtol)
33: !
34: ! End of subroutine FORLINK
35: !
36: return
37: end
39: !============================== Block_Initialization ======== =======72
40: !
41: ! Initializes the common blocks members for turbulence model
42: !
43: !=====================================================================72
44: block data Block_Initialization
45: implicit none
46: #include <petsc/finclude/petscsys.h>
47: PetscScalar vkar,cmu,ce1,ce2
48: PetscScalar aplus1,aplus2,turbinf
49: PetscScalar cb1,sig,cb2,cw1,cw2
50: PetscScalar cw3,cv1,ct1,ct2,ct3,ct4
51: common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
52: common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
53: data vkar,cmu,ce1,ce2/0.41,0.09,1.2,2.0 /
54: data aplus1,aplus2,turbinf/26.0,10.0,0.1 /
55: data cb1,sig,cb2,cw2,cw3/0.1355,0.66667,0.622,0.3,2.0/
56: !
57: ! Comment out old coefficients
58: ! data cv1,ct1,ct2,ct3,ct4/7.1,1.0,2.0,1.1,2.0/
59: !
60: data cv1,ct1,ct2,ct3,ct4/7.1d0,1.0d0,2.0d0,1.2d0,0.5d0/
62: end
64: !================================== INIT =============================72
65: !
66: ! Initializes the flow field
67: !
68: !=====================================================================72
69: subroutine INIT(nnodes, qvec, turbre, amut,nvnode, ivnode,irank)
70: implicit none
71: #include <petsc/finclude/petscsys.h>
72: PetscScalar turbre(1),amut(1)
73: integer ivnode(1),nnodes,nvnode,irank
74: PetscScalar title(20),beta,alpha
75: PetscScalar Re,dt,tot,res0,resc
76: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
77: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
78: PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
79: PetscScalar vkar,cmu,ce1,ce2
80: PetscScalar aplus1,aplus2,turbinf
81: PetscScalar cb1,sig,cb2,cw1,cw2,cw3
82: PetscScalar cv1,ct1,ct2,ct3,ct4
83: PetscScalar cfl1,cfl2
84: integer nsmoth,iflim,itran,nbtran,jupdate, &
85: & nstage,ncyct,iramp,nitfo
86: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
87: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
88: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
89: common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
90: common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
91: common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
92: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
93: & nstage,ncyct,iramp,nitfo
94: PetscScalar pi,conv,rmu,chi,fv1
95: Integer i,k,n
96: #if defined(INTERLACING)
97: PetscScalar qvec(4,nnodes)
98: #define qnod(i,j) qvec(i,j)
99: #else
100: PetscScalar qvec(nnodes,4)
101: #define qnod(i,j) qvec(j,i)
102: #endif
103: !
104: cw1 = cb1/vkar/vkar + (1.0d0 + cb2)/sig
105: if (ivisc.gt.4d0)turbinf = 1.341946d0
106: ! if (ivisc.gt.4)turbinf = 0.5
107: if (ivisc.gt.4.and.itran.eq.1)turbinf = 0.1d0
108: !
109: ! Note that for Spalarts model, I use turbinf as the freestream value of
110: ! the dependent variable just as in the Baldwin-Barth model.
111: ! The constant is set so that in the freestream, nu_t=0.009 (1.341946)
112: !
113: ! print *, "I am in INIT"
114: res0 = 1.0d0
115: resc = 1.0d0
117: gamma = 1.0d0
118: gm1 = 1.0d0
119: gp1 = 1.0d0
120: gm1g = 1.0d0
121: gp1g = 1.0d0
122: ggm1 = 1.0d0
124: pi = 4.0d0*datan(1.0d0)
125: conv = 180.0d0/pi
126: p0 = 1.0d0
127: #if defined(CFL3D_AXIS)
128: u0 = cos(alpha/conv)
129: v0 = 0.0d0
130: w0 = sin(alpha/conv)
131: #else
132: u0 = cos(alpha/conv)
133: v0 = sin(alpha/conv)
134: w0 = 0.0d0
135: #endif
136: do n = 1,nnodes
137: qnod(1,n) = p0
138: qnod(2,n) = u0
139: qnod(3,n) = v0
140: qnod(4,n) = w0
141: enddo
142: if (ivisc.eq.3) then
143: do n = 1,nnodes
144: turbre(n) = 0.0d0
145: amut(n) = 0.0d0
146: enddo
147: endif
148: !
149: ! If viscous, zero out the velocity on the surface
150: !
151: ! print *, "Just Before Viscous"
153: do 9010 i = 1,nvnode
154: !
155: ! Compute the velocity normal to the surface
156: !
157: k = ivnode(i)
158: qnod(2,k) = 0.0d0
159: qnod(3,k) = 0.0d0
160: qnod(4,k) = 0.0d0
161: 9010 continue
162: !
163: ! If turbulent, initialize turbre
164: !
165: if (ivisc.eq.3.or.ivisc.eq.4) then
166: if (irank .eq. 0) write(10,110)
167: if (irank .eq. 0) write(10,120)vkar,cmu,ce1,ce2
168: if (irank .eq. 0) write(10,130)aplus1,aplus2,turbinf
169: do 1010 n = 1,nnodes
170: turbre(n) = turbinf
171: amut(n) = cmu*turbinf
172: 1010 continue
173: end if
175: if (ivisc.eq.5.or.ivisc.eq.6) then
176: if (irank .eq. 0) write(10,110)
177: if (irank .eq. 0) write(10,140)vkar,cb1,sig,cb2
178: if (irank .eq. 0) write(10,150)cw1,cw2,cw3,cv1
179: if (irank .eq. 0) write(10,160)ct1,ct2,ct3,ct4
180: do 1020 n = 1,nnodes
181: turbre(n) = turbinf
182: rmu = 1.0d0
183: chi = turbre(n)/rmu
184: fv1 = chi**3/(chi**3 + cv1**3)
185: amut(n) = fv1*turbinf
186: 1020 continue
187: end if
188: ! print *, "I am out of INIT"
190: return
191: 110 format(1h ,'Parameters for turbulence model')
192: 120 format(1h ,'k=',f10.5,' cmu=',f10.5,' ce1=',f10.5,'ce2=',f10.5)
193: 130 format(1h ,'aplus1',f10.5,' aplus2=',f10.5,' turbinf=',f10.5)
194: 140 format(1h ,'k=',f10.5,' cb1=',f10.5,' sig=',f10.5,'cb2=',f10.5)
195: 150 format(1h ,'cw1=',f10.5,' cw2=',f10.5,' cw3=',f10.5,' cv1=',f10.5)
196: 160 format(1h ,'ct1=',f10.5,' ct2=',f10.5,' ct3=',f10.5,' ct4=',f10.5)
197: !
198: ! End of subroutine INIT
199: !
201: end
203: !================================ READR1 ====================================
204: !
205: ! Reads input parameters
206: !
207: !============================================================================
208: subroutine READR1(ileast, irank)
210: implicit none
211: #include <petsc/finclude/petscsys.h>
212: PetscScalar title(20),beta,alpha
213: PetscScalar Re,dt,tot,res0,resc
214: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
215: PetscScalar cfl1,cfl2
216: integer nsmoth,iflim,itran,nbtran,jupdate, &
217: & nstage,ncyct,iramp,nitfo
218: PetscScalar gtol
219: integer icycle,nsrch,ilu0,ifcn
220: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
221: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
222: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
223: & nstage,ncyct,iramp,nitfo
224: common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
225: integer ileast, irank
226: integer i
228: read(7,10)(title(i),i=1,20)
229: if (irank .eq. 0) write(10,11)(title(i),i=1,20)
231: read(7,10)
233: read(7,24)mseq,ihane,ivisc,ileast,iflim,jupdate,ifcn
234: if (irank .eq. 0) write(10,25)mseq,ihane,ivisc
235: if (irank .eq. 0) write(10,28)ileast,iflim,jupdate
237: read(7,10)
239: read(7,12)beta,alpha,Re
240: if (irank .eq. 0) write(10,13)beta,alpha,Re
242: read(7,10)
244: read(7,26)cfl2,dt,irest,itran,nbtran
245: if (irank .eq. 0) write(10,27)cfl2,dt,irest
247: if (irank .eq. 0) then
248: if (ivisc.eq.5.or.ivisc.eq.6) write(10,123)itran,nbtran
249: endif
251: 10 format(20a4)
252: 11 format(1h ,20a4)
253: 12 format(2f10.5,e14.7)
254: 13 format(1h ,'beta = ',f10.5,' Alpha = ',f10.5,' Re = ',e14.7)
255: 24 format(i10,i10,i10,i10,i10,i10,i10)
256: 25 format(1h ,'mseq = ',i3,' ihane = ',i3,' ivisc=',i3)
257: 26 format(2f10.5,3i10)
258: 27 format(1h ,' cfl2= ',e14.7,' dt= ',f10.5,'irest= ',i5)
259: 28 format(1h ,'ileast= ',i5,' iflim= ',i5,' jupdate= ',i5)
260: 123 format(1h ,'itran = ',i5,' nbtran= ',i5)
262: return
263: end
265: !================================ TECFLO =============================72
266: !
267: ! Writes a formatted file that contains an input file for TECPLOT
268: !
269: !=====================================================================72
270: subroutine TECFLO(nnodes, &
271: & nnbound,nvbound,nfbound, &
272: & nnfacet,nvfacet,nffacet, &
273: & nsnode,nvnode,nfnode, &
274: & title,x,y,z,q, &
275: & nnpts,nntet,nvpts,nvtet,nfpts,nftet, &
276: & f2ntn,f2ntv,f2ntf,isnode,ivnode,ifnode, &
277: & irank)
278: !
279: implicit none
280: #include <petsc/finclude/petscsys.h>
281: PetscScalar pinf
282: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
283: PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
284: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
285: common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
286: integer nnodes,nnbound,nvbound,nfbound, &
287: & nnfacet,nvfacet,nffacet, &
288: & nsnode,nvnode,nfnode,irank
289: integer nntet(1),nnpts(1)
290: integer nvtet(1),nvpts(1)
291: integer nftet(1),nfpts(1)
292: integer f2ntn(nnfacet,4),f2ntv(nvfacet,4),f2ntf(nffacet,4)
293: integer isnode(1),ivnode(1),ifnode(1)
294: PetscScalar x(1),y(1),z(1),q(4,nnodes)
295: integer i,j,i1,i2,i3,ic,isp,ist
297: character c4*4,title*20
299: pinf = p0
300: !
301: ! + write tecplot header
302: !
303: write(13,'(a)') 'TITLE="' // title // '"'
304: write(13,'(a)') 'VARIABLES="X ","Y ","Z ","RHO ",', &
305: & '"U ","V ","W ","P/Pinf"'
306: !
307: ! + do a zone-title, so we can keep track of things
308: ! + start with solid-wall boundary surfaces
309: !
310: isp = 1
311: ist = 1
312: do 10 i=1,nnbound
313: write(c4,"(i4)") i
314: if (i .ge. 0 .and. i .le. 9) ic = 4
315: if (i .ge. 10 .and. i .le. 99) ic = 3
316: if (i .ge. 100 .and. i .le. 999) ic = 2
317: if (i .ge. 1000 .and. i .le. 9999) ic = 1
319: write(13,1000) "nn." // c4(ic:4),nnpts(i),nntet(i)
321: write(13,1010) (x(isnode(j)) ,j=isp,isp+nnpts(i)-1)
322: write(13,1010) (y(isnode(j)) ,j=isp,isp+nnpts(i)-1)
323: write(13,1010) (z(isnode(j)) ,j=isp,isp+nnpts(i)-1)
324: write(13,1010) (1.0 ,j=isp,isp+nnpts(i)-1)
325: write(13,1010) (q(2,isnode(j)) ,j=isp,isp+nnpts(i)-1)
326: write(13,1010) (q(3,isnode(j)) ,j=isp,isp+nnpts(i)-1)
327: write(13,1010) (q(4,isnode(j)) ,j=isp,isp+nnpts(i)-1)
328: write(13,1010) (q(1,isnode(j))/pinf ,j=isp,isp+nnpts(i)-1)
330: do 30 j=ist,ist+nntet(i)-1
331: i1 = f2ntn(j,1) - isp + 1
332: i2 = f2ntn(j,2) - isp + 1
333: i3 = f2ntn(j,3) - isp + 1
335: write(13,1020) i1,i2,i3,i3
336: 30 continue
338: isp = isp + nnpts(i)
339: ist = ist + nntet(i)
341: 10 continue
342: !
343: ! + start with viscous boundary surfaces
344: !
345: isp = 1
346: ist = 1
347: do 70 i=1,nvbound
348: write(c4,"(i4)") i
349: if (i .ge. 0 .and. i .le. 9) ic = 4
350: if (i .ge. 10 .and. i .le. 99) ic = 3
351: if (i .ge. 100 .and. i .le. 999) ic = 2
352: if (i .ge. 1000 .and. i .le. 9999) ic = 1
354: write(13,1000) "nv." // c4(ic:4),nvpts(i),nvtet(i)
356: write(13,1010) (x(ivnode(j)) ,j=isp,isp+nvpts(i)-1)
357: write(13,1010) (y(ivnode(j)) ,j=isp,isp+nvpts(i)-1)
358: write(13,1010) (z(ivnode(j)) ,j=isp,isp+nvpts(i)-1)
359: write(13,1010) (1.0 ,j=isp,isp+nvpts(i)-1)
360: write(13,1010) (q(2,ivnode(j)) ,j=isp,isp+nvpts(i)-1)
361: write(13,1010) (q(3,ivnode(j)) ,j=isp,isp+nvpts(i)-1)
362: write(13,1010) (q(4,ivnode(j)) ,j=isp,isp+nvpts(i)-1)
363: write(13,1010) (q(1,ivnode(j))/pinf ,j=isp,isp+nvpts(i)-1)
365: do 90 j=ist,ist+nvtet(i)-1
366: i1 = f2ntv(j,1) - isp + 1
367: i2 = f2ntv(j,2) - isp + 1
368: i3 = f2ntv(j,3) - isp + 1
370: write(13,1020) i1,i2,i3,i3
371: 90 continue
373: isp = isp + nvpts(i)
374: ist = ist + nvtet(i)
376: 70 continue
377: !
378: ! + do the far-field boundary surfaces
379: !
380: isp = 1
381: ist = 1
382: do 40 i=1,nfbound
383: write(c4,"(i4)") i
384: if (i .ge. 0 .and. i .le. 9) ic = 4
385: if (i .ge. 10 .and. i .le. 99) ic = 3
386: if (i .ge. 100 .and. i .le. 999) ic = 2
387: if (i .ge. 1000 .and. i .le. 9999) ic = 1
389: write(13,1000) "ff." // c4(ic:4),nfpts(i),nftet(i)
391: write(13,1010) (x(ifnode(j)) ,j=isp,isp+nfpts(i)-1)
392: write(13,1010) (y(ifnode(j)) ,j=isp,isp+nfpts(i)-1)
393: write(13,1010) (z(ifnode(j)) ,j=isp,isp+nfpts(i)-1)
394: write(13,1010) (1.0 ,j=isp,isp+nfpts(i)-1)
395: write(13,1010) (q(2,ifnode(j)) ,j=isp,isp+nfpts(i)-1)
396: write(13,1010) (q(3,ifnode(j)) ,j=isp,isp+nfpts(i)-1)
397: write(13,1010) (q(4,ifnode(j)) ,j=isp,isp+nfpts(i)-1)
398: write(13,1010) (q(1,ifnode(j))/pinf ,j=isp,isp+nfpts(i)-1)
400: do 60 j=ist,ist+nftet(i)-1
401: i1 = f2ntf(j,1) - isp + 1
402: i2 = f2ntf(j,2) - isp + 1
403: i3 = f2ntf(j,3) - isp + 1
405: write(13,1020) i1,i2,i3,i3
406: 60 continue
408: isp = isp + nfpts(i)
409: ist = ist + nftet(i)
411: 40 continue
413: !
414: ! End of subroutine TECFLO
415: !
416: 1000 format('ZONE T="',a,'", Z=0, I=',i6,', J=',i6,', F=FEBLOCK')
417: 1010 format(1P10E13.5)
418: 1020 format(4I10)
420: return
421: end
423: !================================ FORCE =============================72
424: !
425: ! Calculates the forces
426: ! Modified - D. K. Kaushik (1/15/97)
427: ! Added new parameters - clift, cdrag, cmom, irank, nvertices
428: !
429: !=====================================================================72
430: subroutine FORCE(nnodes,nedge, &
431: & isnode,ivnode, &
432: & nnfacet,f2ntn,nnbound, &
433: & nvfacet,f2ntv,nvbound, &
434: & evec,qvec,xyz,sface_bit,vface_bit, &
435: & clift, cdrag, cmom,irank,nvertices)
436: implicit none
437: #include <petsc/finclude/petscsys.h>
438: PetscScalar title(20),beta,alpha
439: PetscScalar Re,dt,tot,res0,resc
440: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
441: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
442: PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
443: PetscScalar rms(1001),clw(1001),cdw(1001), &
444: & cmw(1001),xres(1001)
445: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
446: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
447: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
448: common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
449: common/history/rms,clw,cdw,cmw,xres
451: integer nnodes,nedge,nnfacet,nvfacet, &
452: & nnbound,nvbound, &
453: & irank,nvertices
454: integer isnode(1),ivnode(1)
455: integer f2ntn(nnfacet,4)
456: integer f2ntv(nvfacet,4)
457: integer sface_bit(nnfacet), vface_bit(nvfacet)
459: PetscScalar clift, cdrag, cmom
460: PetscScalar clloc,cdloc,cmloc
461: PetscScalar clglo,cdglo,cmglo
462: PetscScalar cl_ghost, cd_ghost, cm_ghost
463: PetscScalar coef_loc(3), coef_glo(3)
464: PetscScalar pi,conv,csa,sna,xr,yr,zr
465: PetscScalar x1,y1,z1,x2,y2,z2,x3,y3
466: PetscScalar z3,xmid,ymid,zmid
467: PetscScalar p1,p2,p3,u1,u2,u3,v1,v2,v3
468: PetscScalar w1,w2,w3,ax,ay,az,bx,by,bz
469: PetscScalar xnorm,ynorm,znorm
470: PetscScalar dcx,dcy,dcz,cp,press
471: integer icount,n,node1,node2,node3,ierr
473: #if defined(INTERLACING)
474: PetscScalar qvec(4,nvertices)
475: PetscScalar xyz(3,nvertices)
476: integer evec(2,nedge)
477: #define qnod(i,j) qvec(i,j)
478: #define x(i) xyz(1,i)
479: #define y(i) xyz(2,i)
480: #define z(i) xyz(3,i)
481: #define eptr(j,i) evec(i,j)
482: #else
483: PetscScalar qvec(nvertices,4)
484: PetscScalar xyz(nvertices,3)
485: integer evec(nedge,2)
486: #define qnod(i,j) qvec(j,i)
487: #define x(i) xyz(i,1)
488: #define y(i) xyz(i,2)
489: #define z(i) xyz(i,3)
490: #define eptr(i,j) evec(i,j)
491: #endif
492: !
493: pi = 4.d0*atan(1.d0)
494: conv = 180.d0/pi
495: csa=cos(alpha/conv)
496: sna=sin(alpha/conv)
497: !
498: ! initialize forces to zero
499: !
500: clw(ntt) = 0.0d0
501: cdw(ntt) = 0.0d0
502: cmw(ntt) = 0.0d0
503: cl_ghost = 0.0d0
504: cd_ghost = 0.0d0
505: cm_ghost = 0.0d0
506: clloc = 0.0d0
507: cdloc = 0.0d0
508: cmloc = 0.0d0
509: clglo = 0.0d0
510: cdglo = 0.0d0
511: cmglo = 0.0d0
512: coef_glo(1) = 0.0d0
513: coef_glo(2) = 0.0d0
514: coef_glo(3) = 0.0d0
516: xr = 0.25d0
517: yr = 0.0d0
518: zr = 0.0d0
520: do 30 n = 1, nnfacet
521: node1 = isnode(f2ntn(n,1))
522: node2 = isnode(f2ntn(n,2))
523: node3 = isnode(f2ntn(n,3))
525: x1 = x(node1)
526: y1 = y(node1)
527: z1 = z(node1)
529: x2 = x(node2)
530: y2 = y(node2)
531: z2 = z(node2)
533: x3 = x(node3)
534: y3 = y(node3)
535: z3 = z(node3)
537: ax = x2 - x1
538: ay = y2 - y1
539: az = z2 - z1
541: bx = x3 - x1
542: by = y3 - y1
543: bz = z3 - z1
544: !
545: ! norm points outward, away from grid interior.
546: ! norm magnitude is area of surface triangle.
547: !
548: xnorm =-0.5d0*(ay*bz - az*by)
549: ynorm = 0.5d0*(ax*bz - az*bx)
550: znorm =-0.5d0*(ax*by - ay*bx)
552: p1 = qnod(1,node1)
553: u1 = qnod(2,node1)
554: v1 = qnod(3,node1)
555: w1 = qnod(4,node1)
557: p2 = qnod(1,node2)
558: u2 = qnod(2,node2)
559: v2 = qnod(3,node2)
560: w2 = qnod(4,node2)
562: p3 = qnod(1,node3)
563: u3 = qnod(2,node3)
564: v3 = qnod(3,node3)
565: w3 = qnod(4,node3)
567: press = (p1 + p2 + p3)/3.0d0
568: cp = 2.d0*(press - 1.0d0)
569: !
570: dcx = cp*xnorm
571: dcy = cp*ynorm
572: dcz = cp*znorm
574: xmid = (x1 + x2 + x3)
575: ymid = (y1 + y2 + y3)
576: zmid = (z1 + z2 + z3)
578: if (sface_bit(n) .eq. 1) then
579: #if defined(CFL3D_AXIS)
581: clloc = clloc - dcx*sna + dcz*csa
582: cdloc = cdloc + dcx*csa + dcz*sna
583: cmloc = cmloc &
584: & + (xmid - xr)*dcy - (ymid - yr)*dcx
586: #else
587: clloc = clloc - dcx*sna + dcy*csa
588: cdloc = cdloc + dcx*csa + dcy*sna
589: cmloc = cmloc &
590: & + (xmid - xr)*dcy - (ymid - yr)*dcx
591: #endif
592: endif
594: 30 continue
595: !
596: ! Viscous boundary
597: !
598: do 60 n = 1, nvfacet
599: node1 = ivnode(f2ntv(n,1))
600: node2 = ivnode(f2ntv(n,2))
601: node3 = ivnode(f2ntv(n,3))
603: x1 = x(node1)
604: y1 = y(node1)
605: z1 = z(node1)
607: x2 = x(node2)
608: y2 = y(node2)
609: z2 = z(node2)
611: x3 = x(node3)
612: y3 = y(node3)
613: z3 = z(node3)
615: ax = x2 - x1
616: ay = y2 - y1
617: az = z2 - z1
619: bx = x3 - x1
620: by = y3 - y1
621: bz = z3 - z1
622: !
623: ! norm points outward, away from grid interior.
624: ! norm magnitude is area of surface triangle.
625: !
626: xnorm =-0.5d0*(ay*bz - az*by)
627: ynorm = 0.5d0*(ax*bz - az*bx)
628: znorm =-0.5d0*(ax*by - ay*bx)
630: p1 = qnod(1,node1)
631: u1 = qnod(2,node1)
632: v1 = qnod(3,node1)
633: w1 = qnod(4,node1)
635: p2 = qnod(1,node2)
636: u2 = qnod(2,node2)
637: v2 = qnod(3,node2)
638: w2 = qnod(4,node2)
640: p3 = qnod(1,node3)
641: u3 = qnod(2,node3)
642: v3 = qnod(3,node3)
643: w3 = qnod(4,node3)
645: press = (p1 + p2 + p3)/3.0d0
646: cp = 2.d0*(press -1.0d0)
648: dcx = cp*xnorm
649: dcy = cp*ynorm
650: dcz = cp*znorm
652: xmid = (x1 + x2 + x3)
653: ymid = (y1 + y2 + y3)
654: zmid = (z1 + z2 + z3)
656: if (vface_bit(n) .eq. 1) then
657: clloc = clloc - dcx*sna + dcy*csa
658: cdloc = cdloc + dcx*csa + dcy*sna
659: cmloc = cmloc &
660: & + (xmid - xr)*dcy - (ymid - yr)*dcx
661: endif
663: 60 continue
664: !
665: icount = 3
666: coef_loc(1) = clloc
667: coef_loc(2) = cdloc
668: coef_loc(3) = cmloc
669: call MPI_ALLREDUCE(coef_loc,coef_glo,icount, &
670: & MPIU_SCALAR, MPIU_SUM, &
671: & MPI_COMM_WORLD,ierr)
672: ! call MPI_ALLREDUCE(cdloc,cdglo,icount,
673: ! & MPIU_REAL_PRECISION, MPI_SUM,
674: ! & MPI_COMM_WORLD,ierr)
675: ! call MPI_ALLREDUCE(cmloc,cmglo,icount,
676: ! & MPIU_REAL_PRECISION, MPI_SUM,
677: ! & MPI_COMM_WORLD,ierr)
678: !
679: clw(ntt) = coef_glo(1)
680: cdw(ntt) = coef_glo(2)
681: cmw(ntt) = coef_glo(3)
682: clift = clw(ntt)
683: cdrag = cdw(ntt)
684: cmom = cmw(ntt)
685: !
686: ! Update the timing information - Added by D. K. Kaushik (1/17/97)
687: xres(ntt) = tot
688: !
689: ! End of subroutine FORCE
690: !
691: return
692: end
694: !================================ DELTAT2 ============================72
695: !
696: ! Calculate a time step for each cell
697: ! Note that this routine assumes conservative variables
698: !
699: !=====================================================================72
700: subroutine DELTAT2(nnodes,nedge,qvec,cdt, &
701: & xyz,vol,xyzn,evec, &
702: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
703: & nsnode,nvnode,nfnode,isnode,ivnode,ifnode, &
704: & LocalTS,irank,nvertices)
705: !
706: implicit none
707: #include <petsc/finclude/petscsys.h>
708: PetscScalar title(20),beta,alpha
709: PetscScalar Re,dt,tot,res0,resc
710: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
711: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
712: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
713: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
714: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
715: integer nnodes,nedge,nsnode,nvnode,nfnode,LocalTS,irank,nvertices
716: PetscScalar vol(1)
717: PetscScalar sxn(1),syn(1),szn(1)
718: PetscScalar vxn(1),vyn(1),vzn(1)
719: PetscScalar fxn(1),fyn(1),fzn(1)
720: PetscScalar cdt(1)
721: integer isnode(1),ivnode(1),ifnode(1)
722: integer n,i,node1,node2,inode
723: integer ierr
724: PetscLogDouble flops
725: PetscScalar xnorm,ynorm,znorm,area
726: PetscScalar u1,v1,w1,u2,v2,w2,u,v,c,w
727: PetscScalar ubar,Vn,term
728: #if defined(INTERLACING)
729: PetscScalar qvec(4,nvertices)
730: PetscScalar xyz(3,nvertices)
731: PetscScalar xyzn(4,nedge)
732: integer evec(2,nedge)
733: #define qnod(i,j) qvec(i,j)
734: #define x(i) xyz(1,i)
735: #define y(i) xyz(2,i)
736: #define z(i) xyz(3,i)
737: #define xn(i) xyzn(1,i)
738: #define yn(i) xyzn(2,i)
739: #define zn(i) xyzn(3,i)
740: #define rl(i) xyzn(4,i)
741: #define eptr(j,i) evec(i,j)
742: #else
743: PetscScalar qvec(nvertices,4)
744: PetscScalar xyz(nvertices,3)
745: PetscScalar xyzn(nedge,4)
746: integer evec(nedge,2)
747: #define qnod(i,j) qvec(j,i)
748: #define x(i) xyz(i,1)
749: #define y(i) xyz(i,2)
750: #define z(i) xyz(i,3)
751: #define xn(i) xyzn(i,1)
752: #define yn(i) xyzn(i,2)
753: #define zn(i) xyzn(i,3)
754: #define rl(i) xyzn(i,4)
755: #define eptr(i,j) evec(i,j)
756: #endif
757: !
758: ! If local time steping, loop over faces
759: ! and calculate time step as cdt = V/(sum(|u.n| +c.area)
760: ! This is time step for cfl=1. We will multiply by cfl number later
761: !
762: ! First loop over nodes and zero out cdt
763: !
764: flops = 0.0
765: if (LocalTS .gt. 0) then
766: do 1000 i = 1,nvertices
767: cdt(i) = 0.0d0
768: 1000 continue
769: !
770: do 1020 n = 1, nedge
771: node1 = eptr(n,1)
772: node2 = eptr(n,2)
773: !
774: ! Get normal to face
775: !
776: xnorm = xn(n)
777: ynorm = yn(n)
778: znorm = zn(n)
779: area = rl(n)
780: !
781: xnorm = xnorm*area
782: ynorm = ynorm*area
783: znorm = znorm*area
784: !
785: !/*
786: ! xnorm = xnormal * area of face
787: ! ynorm = ynormal * area of face
788: ! znorm = znormal * area of face
789: !*/
790: u1 = qnod(2,node1)
791: v1 = qnod(3,node1)
792: w1 = qnod(4,node1)
794: u2 = qnod(2,node2)
795: v2 = qnod(3,node2)
796: w2 = qnod(4,node2)
797: !
798: ! Get average values on face
799: !
800: u = 0.5d0*(u1 + u2)
801: v = 0.5d0*(v1 + v2)
802: w = 0.5d0*(w1 + w2)
803: ubar = xn(n)*u + yn(n)*v + zn(n)*w
804: c = sqrt(ubar*ubar + beta)
806: term = abs(u*xnorm + v*ynorm + w*znorm) + c*area
807: cdt(node1) = cdt(node1) + term
808: cdt(node2) = cdt(node2) + term
809: !
810: 1020 continue
811: flops = flops + 27.0*nedge
812: !
813: ! Now loop over boundaries and close the contours
814: !
815: do 1030 i = 1,nsnode
816: inode = isnode(i)
817: !
818: ! Get the normal
819: !
820: xnorm = sxn(i)
821: ynorm = syn(i)
822: znorm = szn(i)
823: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
825: u = qnod(2,inode)
826: v = qnod(3,inode)
827: w = qnod(4,inode)
828: ubar= (u*xnorm + v*ynorm + w*znorm)/area
829: c = sqrt(ubar*ubar + beta)
830: !
831: Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
832: cdt(inode) = cdt(inode) + Vn
833: !
834: 1030 continue
835: flops = flops + 24.0*nsnode
837: !
838: ! Now viscous faces
839: !
840: !
841: do 1040 i = 1,nvnode
842: inode = ivnode(i)
843: !
844: ! Get the normal
845: !
846: xnorm = vxn(i)
847: ynorm = vyn(i)
848: znorm = vzn(i)
849: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
850: !
851: ! Because velocities are zero, c is just sqrt(beta)
852: !
853: c = sqrt(beta)
854: Vn = c*area
856: cdt(inode) = cdt(inode) + Vn
857: !
858: 1040 continue
859: flops = flops + 9.0*nvnode
860: !
861: ! Now far field
862: !
863: do 1050 i = 1,nfnode
864: inode = ifnode(i)
865: !
866: ! Get the normal
867: !
868: xnorm = fxn(i)
869: ynorm = fyn(i)
870: znorm = fzn(i)
871: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
873: u = qnod(2,inode)
874: v = qnod(3,inode)
875: w = qnod(4,inode)
876: ubar= (u*xnorm + v*ynorm + w*znorm)/area
877: c = sqrt(ubar*ubar + beta)
879: Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
880: cdt(inode) = cdt(inode) + Vn
881: 1050 continue
882: flops = flops + 24.0*nfnode
884: do 1060 n = 1,nvertices
885: cdt(n) = vol(n)/cdt(n)
886: 1060 continue
887: flops = flops + nvertices
888: else
889: !
890: ! If not doing local time stepping just set cdt=1
891: !
892: do 1070 n = 1,nvertices
893: cdt(n) = 1.0d0
894: 1070 continue
895: end if
896: call PetscLogFlops(flops,ierr)
897: !
898: ! End of subroutine DELTAT2
899: !
900: return
901: end
903: !================================= FLUX ======================================
904: !
905: ! Calculates the fluxes on the face and performs the flux balance
906: !
907: !=============================================================================
908: #if defined(_OPENMP)
909: #if defined(HAVE_EDGE_COLORING)
910: subroutine FLUX(nnodes, ncell, nedge, &
911: & nsface, nvface, nfface, isface, ivface, ifface, &
912: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
913: & nnfacet,f2ntn,nnbound, &
914: & nvfacet,f2ntv,nvbound, &
915: & nffacet,f2ntf,nfbound, &
916: & grad,evec,qvec,xyz,resvec,xyzn, &
917: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
918: & max_threads, &
919: & ncolor,ncount, &
920: & irank,nvertices)
921: #elif defined(HAVE_REDUNDANT_WORK)
922: subroutine FLUX(nnodes, ncell, nedge, &
923: & nsface, nvface, nfface, isface, ivface, ifface, &
924: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
925: & nnfacet,f2ntn,nnbound, &
926: & nvfacet,f2ntv,nvbound, &
927: & nffacet,f2ntf,nfbound, &
928: & grad,evec,qvec,xyz,resvec,xyzn, &
929: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
930: & max_threads, &
931: & resd, &
932: & irank,nvertices)
933: #else
934: subroutine FLUX(nnodes, ncell, nedge, &
935: & nsface, nvface, nfface, isface, ivface, ifface, &
936: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
937: & nnfacet,f2ntn,nnbound, &
938: & nvfacet,f2ntv,nvbound, &
939: & nffacet,f2ntf,nfbound, &
940: & grad,evec,qvec,xyz,resvec,xyzn, &
941: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
942: & max_threads, &
943: & nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr, &
944: & irank,nvertices)
945: #endif
946: #else
947: subroutine FLUX(nnodes, ncell, nedge, &
948: & nsface, nvface, nfface, isface, ivface, ifface, &
949: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
950: & nnfacet,f2ntn,nnbound, &
951: & nvfacet,f2ntv,nvbound, &
952: & nffacet,f2ntf,nfbound, &
953: & grad,evec,qvec,xyz,resvec,xyzn, &
954: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
955: & irank,nvertices)
956: #endif
957: implicit none
958: #include <petsc/finclude/petscsys.h>
959: PetscScalar title(20),beta,alpha
960: PetscScalar Re,dt,tot,res0,resc
961: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
962: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
963: PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
964: PetscScalar rms(1001),clw(1001),cdw(1001)
965: PetscScalar cmw(1001),xres(1001)
966: PetscScalar cfl1,cfl2
967: integer nsmoth,iflim,itran,nbtran,jupdate, &
968: & nstage,ncyct,iramp,nitfo
969: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
970: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
971: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
972: common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
973: common/history/rms,clw,cdw,cmw,xres
974: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
975: & nstage,ncyct,iramp,nitfo
976: integer nnodes, ncell, nedge, &
977: & nsface, nvface, nfface, &
978: & nsnode, nvnode, nfnode, &
979: & nnfacet,nvfacet,nffacet, &
980: & nnbound,nvbound,nfbound, &
981: & irank,nvertices
982: PetscScalar sxn(1),syn(1),szn(1)
983: PetscScalar vxn(1),vyn(1),vzn(1)
984: PetscScalar fxn(1),fyn(1),fzn(1)
985: PetscScalar phi(nvertices,4)
986: integer isface(1),ivface(1),ifface(1)
987: integer isnode(1),ivnode(1),ifnode(1)
988: integer f2ntn(nnfacet,4)
989: integer f2ntv(nvfacet,4)
990: integer f2ntf(nffacet,4)
991: integer i,n,node1,node2,node3,inode
992: integer ierr,first_edge,last_edge
993: PetscLogDouble flops
994: PetscScalar xmean,ymean,zmean,xnorm,ynorm
995: PetscScalar znorm,rlen,area,dot,size
996: PetscScalar X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3
997: PetscScalar p1,p2,p3
998: PetscScalar rx,ry,rz,pL,uL,wL,vL
999: PetscScalar pR,uR,vR,wR
1000: PetscScalar ubarL,ubarR,ubar,p,u,v,w,c
1001: PetscScalar c2,c20
1002: PetscScalar ubar0,pi,ui,vi,wi,unorm
1003: PetscScalar phi1,phi2,phi3,phi4,phi5
1004: PetscScalar phi6,phi7,phi8,phi9
1005: PetscScalar eig1,eig2,eig3,eig4
1006: PetscScalar dp,du,dv,dw
1007: PetscScalar ti11,ti21,ti31,ti41
1008: PetscScalar dv1,dv2,dv3,dv4
1009: PetscScalar r11,r21,r31,r41
1010: PetscScalar r12,r22,r32,r42
1011: PetscScalar r13,r23,r33,r43
1012: PetscScalar r14,r24,r34,r44
1013: PetscScalar t1,t2,t3,t4
1014: PetscScalar t11,t21,t31,t41
1015: PetscScalar t12,t22,t32,t42
1016: PetscScalar t13,t23,t33,t43
1017: PetscScalar t14,t24,t34,t44
1018: PetscScalar fluxp1,fluxp2,fluxp3,fluxp4
1019: PetscScalar fluxm1,fluxm2,fluxm3,fluxm4
1020: PetscScalar res1,res2,res3,res4
1021: PetscScalar rhs1,rhs2,rhs3,rhs4
1022: PetscScalar c68,c18,second
1023: PetscScalar ax,ay,az,bx,by,bz
1024: PetscScalar pa,pb,pc,ub,vb,wb
1025: #if defined(INTERLACING)
1026: PetscScalar qvec(4,nvertices)
1027: PetscScalar grad(3,4,nvertices)
1028: PetscScalar resvec(4,nnodes)
1029: PetscScalar xyz(3,nvertices)
1030: PetscScalar xyzn(4,nedge)
1031: integer evec(2,nedge)
1032: #define qnod(i,j) qvec(i,j)
1033: #define res(i,j) resvec(i,j)
1034: #define gradx(i,j) grad(1,i,j)
1035: #define grady(i,j) grad(2,i,j)
1036: #define gradz(i,j) grad(3,i,j)
1037: #define x(i) xyz(1,i)
1038: #define y(i) xyz(2,i)
1039: #define z(i) xyz(3,i)
1040: #define xn(i) xyzn(1,i)
1041: #define yn(i) xyzn(2,i)
1042: #define zn(i) xyzn(3,i)
1043: #define ra(i) xyzn(4,i)
1044: #define eptr(j,i) evec(i,j)
1045: #else
1046: PetscScalar qvec(nvertices,4)
1047: PetscScalar grad(nvertices,4,3)
1048: PetscScalar resvec(nnodes,4)
1049: PetscScalar xyz(nvertices,3)
1050: PetscScalar xyzn(nedge,4)
1051: integer evec(nedge,2)
1052: #define qnod(i,j) qvec(j,i)
1053: #define res(i,j) resvec(j,i)
1054: #define gradx(i,j) grad(j,i,1)
1055: #define grady(i,j) grad(j,i,2)
1056: #define gradz(i,j) grad(j,i,3)
1057: #define x(i) xyz(i,1)
1058: #define y(i) xyz(i,2)
1059: #define z(i) xyz(i,3)
1060: #define xn(i) xyzn(i,1)
1061: #define yn(i) xyzn(i,2)
1062: #define zn(i) xyzn(i,3)
1063: #define ra(i) xyzn(i,4)
1064: #define eptr(i,j) evec(i,j)
1065: #endif
1066: #if defined(_OPENMP)
1067: integer my_thread_id,omp_get_thread_num, &
1068: & max_threads,omp_get_num_threads, &
1069: & first_node,last_node,chunk
1070: #if defined(HAVE_EDGE_COLORING)
1071: integer ncolor,ncount(1)
1072: #elif defined(HAVE_REDUNDANT_WORK)
1073: PetscScalar resd(4,nnodes)
1074: #undef res
1075: #define res(a,b) resd(a,b)
1076: #else
1077: integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads), &
1078: & edge_thr(2,nedgeAllThr)
1079: PetscScalar xyzn_thr(4,nedgeAllThr)
1080: #undef eptr
1081: #define eptr(a,b) edge_thr(b,a)
1082: #undef xn
1083: #undef yn
1084: #undef zn
1085: #undef ra
1086: #define xn(i) xyzn_thr(1,i)
1087: #define yn(i) xyzn_thr(2,i)
1088: #define zn(i) xyzn_thr(3,i)
1089: #define ra(i) xyzn_thr(4,i)
1090: #endif
1091: #endif
1092: ! logging variables
1093: integer flux_event, fun3d_class, flag
1094: character * 16 flux_label
1095: data flag/-1/,flux_label/'FluxEvaluation '/
1096: save flux_event, fun3d_class, flag, flux_label
1097: if (flag .eq. -1) then
1098: call PetscClassIdRegister('PetscFun3d',fun3d_class,ierr)
1099: call PetscLogEventRegister(flux_label, &
1100: & fun3d_class,flux_event,ierr)
1101: flag = 1
1102: endif
1103: call PetscLogEventBegin(flux_event,ierr)
1104: flops = 0.0
1105: second = 1.0d0
1106: ! if (ntt.le.nitfo)second = 0.0
1107: ! print *, 'Second is' , second
1108: #if defined(_OPENMP)
1109: #if defined(HAVE_EDGE_COLORING)
1110: first_edge = 1
1111: do i = 1,ncolor
1112: last_edge = first_edge + ncount(i) - 1
1113: #elif defined(HAVE_REDUNDANT_WORK)
1114: do i = 1,nnodes
1115: do j = 1,4
1116: resd(j,i) = 0.0
1117: enddo
1118: enddo
1119: first_edge = 1
1120: last_edge = nedge
1121: #endif
1122: #else
1123: first_edge = 1
1124: last_edge = nedge
1125: #endif
1126: !
1127: ! Loop over all the faces and calculate flux
1128: !
1129: !$omp parallel default(shared) private(n,node1,node2,xmean,ymean, &
1130: !$omp& zmean,xnorm,ynorm,znorm,rlen,dot,x1,y1,z1,size,x2,y2,z2,rx, &
1131: !$omp& ry,rz,pL,uL,vL,wL,ubarL,pR,uR,vR,wR,ubarR,p,u,v,w,ubar,phi1, &
1132: !$omp& phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,c2,c,eig1,eig2,eig3, &
1133: !$omp& eig4,dp,du,dv,dw,ti11,ti21,ti31,ti41,dv1,dv2,dv3,dv4,r11,r21, &
1134: !$omp& r31,r41,r12,r22,r32,r42,r13,r23,r33,r43,r14,r24,r34,r44,t1,t2, &
1135: !$omp& t3,t4,fluxp1,fluxp2,fluxp3,fluxp4,fluxm1,fluxm2,fluxm3,fluxm4, &
1136: !$omp& res1,res2,res3,res4,c68,c18,inode) &
1137: !$omp& private(x3,y3,z3,p1,p2,p3,ax,ay,az,bx,by,bz,pa,pb,pc,area,c20, &
1138: !$omp& c0,ubar0,t11,t21,t31,t41,t12,t22,t32,t42,t13,t23,t33,t43,t14, &
1139: !$omp& t24,t34,t44,pi,ui,vi,wi,unorm,rhs1,rhs2,rhs3,rhs4,ub,vb,wb,node3) &
1140: !$omp& private(my_thread_id,chunk) &
1141: !$omp& reduction(+:flops) &
1142: #if defined(_OPENMP)
1143: #if defined(HAVE_EDGE_COLORING)
1144: #elif defined(HAVE_REDUNDANT_WORK)
1145: !$omp& firstprivate(resd)
1146: #else
1147: !$omp& private(first_edge,last_edge)
1148: #endif
1149: #endif
1150: !!$omp1 shared(evec,qvec,resvec,xyz,xyzn,grad,nedge,nnodes,second,flops,
1151: !!$omp2 title,beta,alpha,Re,dt,tot,res0,resc,ntt,mseq,ivisc,irest,icyc,
1152: !!$omp3 ihane,ntturb,gamma,gm1,gp1,gm1g,gp1g,ggm1,p0,rho0,c0,u0,v0,w0,
1153: !!$omp4 et0,h0,pt0,rms,clw,cdw,cmw,xres,cfl1,cfl2,nsmoth,iflim,
1154: !!$omp5 itran,nbtran,jupdate,nstage,ncyct,iramp,nitfo)
1155: !!$ my_thread_id = omp_get_thread_num()
1156: !!$omp critical
1157: !!$ print *, 'My thread id on processor ',irank,' is ',my_thread_id
1158: !!$omp end critical
1159: #if defined(_OPENMP)
1160: #if defined(HAVE_EDGE_COLORING)
1161: !$omp do
1162: #elif defined(HAVE_REDUNDANT_WORK)
1163: !$omp do
1164: #else
1165: ! tot_threads = omp_get_num_threads()
1166: my_thread_id = omp_get_thread_num()
1167: ! chunk = nnodes/tot_threads
1168: ! first_node = my_thread_id*chunk+1
1169: ! if (my_thread_id .eq. (size-1)) then
1170: ! last_node = nnodes
1171: ! else
1172: ! last_node = (my_thread_id+1)*chunk
1173: ! endif
1174: first_edge = nedge_thr(my_thread_id)
1175: last_edge = nedge_thr(my_thread_id+1)-1
1177: ! print *, 'On thread ', my_thread_id,', first_edge = ',first_edge, &
1178: ! & ', last_edge = ',last_edge
1180: #endif
1181: #else
1182: #endif
1183: do n = first_edge, last_edge
1184: node1 = eptr(n,1)
1185: node2 = eptr(n,2)
1186: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1187: if ((part_thr(node1) .eq. my_thread_id) &
1188: & .or.(part_thr(node2) .eq. my_thread_id)) then
1189: #endif
1190: !
1191: ! Calculate unit normal to face and length of face
1192: !
1193: xmean = .5d0*(x(node1) + x(node2))
1194: ymean = .5d0*(y(node1) + y(node2))
1195: zmean = .5d0*(z(node1) + z(node2))
1196: xnorm = xn(n)
1197: ynorm = yn(n)
1198: znorm = zn(n)
1199: rlen = ra(n)
1200: !
1201: ! Now lets get our other 2 vectors
1202: ! For first vector, use {1,0,0} and subtract off the component
1203: ! in the direction of the face normal. If the inner product of
1204: ! {1,0,0} is close to unity, use {0,1,0}
1205: !
1206: dot = xnorm
1207: if (abs(dot).lt.0.95d0)then
1208: X1 = 1. - dot*xnorm
1209: Y1 = - dot*ynorm
1210: Z1 = - dot*znorm
1211: else
1212: dot = ynorm
1213: X1 = - dot*xnorm
1214: Y1 = 1.0d0 - dot*ynorm
1215: Z1 = - dot*znorm
1216: end if
1217: !
1218: ! Normalize the first vector
1219: !
1220: size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
1221: X1 = X1/size
1222: Y1 = Y1/size
1223: Z1 = Z1/size
1224: !
1225: ! Take cross-product of normal and V1 to get V2
1226: !
1227: X2 = ynorm*Z1 - znorm*Y1
1228: Y2 = znorm*X1 - xnorm*Z1
1229: Z2 = xnorm*Y1 - ynorm*X1
1230: !
1231: ! Get variables on "left" and "right" side of face
1232: !
1233: rx = second*(xmean - x(node1))
1234: ry = second*(ymean - y(node1))
1235: rz = second*(zmean - z(node1))
1236: pL = qnod(1,node1) + gradx(1,node1)*rx &
1237: & + grady(1,node1)*ry &
1238: & + gradz(1,node1)*rz
1239: uL = qnod(2,node1) + gradx(2,node1)*rx &
1240: & + grady(2,node1)*ry &
1241: & + gradz(2,node1)*rz
1242: vL = qnod(3,node1) + gradx(3,node1)*rx &
1243: & + grady(3,node1)*ry &
1244: & + gradz(3,node1)*rz
1245: wL = qnod(4,node1) + gradx(4,node1)*rx &
1246: & + grady(4,node1)*ry &
1247: & + gradz(4,node1)*rz
1249: ubarL = xnorm*uL + ynorm*vL + znorm*wL
1250: ! c2L = ubarl*ubarL + beta
1251: ! cL = sqrt(c2L)
1252: !
1253: rx = second*(xmean - x(node2))
1254: ry = second*(ymean - y(node2))
1255: rz = second*(zmean - z(node2))
1256: pR = qnod(1,node2) + gradx(1,node2)*rx &
1257: & + grady(1,node2)*ry &
1258: & + gradz(1,node2)*rz
1259: uR = qnod(2,node2) + gradx(2,node2)*rx &
1260: & + grady(2,node2)*ry &
1261: & + gradz(2,node2)*rz
1262: vR = qnod(3,node2) + gradx(3,node2)*rx &
1263: & + grady(3,node2)*ry &
1264: & + gradz(3,node2)*rz
1265: wR = qnod(4,node2) + gradx(4,node2)*rx &
1266: & + grady(4,node2)*ry &
1267: & + gradz(4,node2)*rz
1268: ubarR = xnorm*uR + ynorm*vR + znorm*wR
1269: ! c2R = ubarR*ubarR + beta
1270: ! cR = sqrt(c2R)
1271: !
1272: ! Compute averages
1273: !
1274: p = .5d0*(pL + pR)
1275: u = .5d0*(uL + uR)
1276: v = .5d0*(vL + vR)
1277: w = .5d0*(wL + wR)
1279: ubar = xnorm*u + ynorm*v + znorm*w
1280: phi1 = xnorm*beta + u*ubar
1281: phi2 = ynorm*beta + v*ubar
1282: phi3 = znorm*beta + w*ubar
1283: phi4 = Y2*phi3 - Z2*phi2
1284: phi5 = Z2*phi1 - X2*phi3
1285: phi6 = X2*phi2 - Y2*phi1
1286: phi7 = Z1*phi2 - Y1*phi3
1287: phi8 = X1*phi3 - Z1*phi1
1288: phi9 = Y1*phi1 - X1*phi2
1289: c2 = ubar*ubar + beta
1290: c = sqrt(c2)
1291: !
1292: ! Now compute eigenvalues, eigenvectors, and strengths
1293: !
1294: eig1 = abs(ubar)
1295: eig2 = abs(ubar)
1296: eig3 = abs(ubar + c)
1297: eig4 = abs(ubar - c)
1298: !
1299: dp = pr - pl
1300: du = ur - ul
1301: dv = vr - vl
1302: dw = wr - wl
1303: !
1304: ! Components of T(inverse) (I will divide by c2 later)
1305: !
1306: ti11 = -(u*phi4 + v*phi5 + w*phi6)/beta
1307: ti21 = -(u*phi7 + v*phi8 + w*phi9)/beta
1308: ti31 = .5d0*(c-ubar)/beta
1309: ti41 = -.5d0*(c+ubar)/beta
1310: !
1311: ! jumps (T(inverse)*dq)
1312: dv1 = (ti11*dp + phi4*du + phi5*dv + phi6*dw)/c2
1313: dv2 = (ti21*dp + phi7*du + phi8*dv + phi9*dw)/c2
1314: dv3 = .5d0*(2.d0*ti31*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
1315: dv4 = .5d0*(2.d0*ti41*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
1316: !
1317: ! Now get elements of T (call it r for now)
1318: !
1319: r11 = 0.0d0
1320: r21 = X1
1321: r31 = Y1
1322: r41 = Z1
1324: r12 = 0.0d0
1325: r22 = X2
1326: r32 = Y2
1327: r42 = Z2
1329: r13 = c*beta
1330: r23 = xnorm*beta + u*(ubar + c)
1331: r33 = ynorm*beta + v*(ubar + c)
1332: r43 = znorm*beta + w*(ubar + c)
1333: !
1334: r14 = -c*beta
1335: r24 = xnorm*beta + u*(ubar - c)
1336: r34 = ynorm*beta + v*(ubar - c)
1337: r44 = znorm*beta + w*(ubar - c)
1338: !
1339: ! Calculate T* |lambda| *T(inverse)
1340: !
1341: t1 = eig1*r11*dv1 + eig2*r12*dv2 + eig3*r13*dv3 + eig4*r14*dv4
1342: t2 = eig1*r21*dv1 + eig2*r22*dv2 + eig3*r23*dv3 + eig4*r24*dv4
1343: t3 = eig1*r31*dv1 + eig2*r32*dv2 + eig3*r33*dv3 + eig4*r34*dv4
1344: t4 = eig1*r41*dv1 + eig2*r42*dv2 + eig3*r43*dv3 + eig4*r44*dv4
1345: !
1346: ! Modify to calculate .5(fl +fr) from nodes
1347: ! instead of extrapolated ones
1348: !
1349: ! pL = qnod(1,node1)
1350: ! uL = qnod(2,node1)
1351: ! vL = qnod(3,node1)
1352: ! wL = qnod(4,node1)
1353: ! ubarL = xnorm*uL + ynorm*vL + znorm*wL
1354: !
1355: fluxp1 = rlen*beta*ubarL
1356: fluxp2 = rlen*(uL*ubarL + xnorm*pL)
1357: fluxp3 = rlen*(vL*ubarL + ynorm*pL)
1358: fluxp4 = rlen*(wL*ubarL + znorm*pL)
1359: !
1360: ! Now the right side
1361: !
1362: ! pR = qnod(1,node2)
1363: ! uR = qnod(2,node2)
1364: ! vR = qnod(3,node2)
1365: ! wR = qnod(4,node2)
1366: ! ubarR = xnorm*uR + ynorm*vR + znorm*wR
1367: !
1368: fluxm1 = rlen*beta*ubarR
1369: fluxm2 = rlen*(uR*ubarR + xnorm*pR)
1370: fluxm3 = rlen*(vR*ubarR + ynorm*pR)
1371: fluxm4 = rlen*(wR*ubarR + znorm*pR)
1372: !
1373: res1 = 0.5d0*(fluxp1 + fluxm1 - rlen*t1)
1374: res2 = 0.5d0*(fluxp2 + fluxm2 - rlen*t2)
1375: res3 = 0.5d0*(fluxp3 + fluxm3 - rlen*t3)
1376: res4 = 0.5d0*(fluxp4 + fluxm4 - rlen*t4)
1377: !
1378: flops = flops + 318.0
1380: !
1381: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1382: if ((part_thr(node1).eq.my_thread_id) &
1383: & .and.(node1.le.nnodes)) then
1384: #else
1385: if (node1.le.nnodes) then
1386: #endif
1387: res(1,node1) = res(1,node1) + res1
1388: res(2,node1) = res(2,node1) + res2
1389: res(3,node1) = res(3,node1) + res3
1390: res(4,node1) = res(4,node1) + res4
1391: flops = flops + 4.0
1392: endif
1393: !
1394: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
1395: if ((part_thr(node2).eq.my_thread_id) &
1396: & .and.(node2.le.nnodes)) then
1397: #else
1398: if (node2.le.nnodes) then
1399: #endif
1400: res(1,node2) = res(1,node2) - res1
1401: res(2,node2) = res(2,node2) - res2
1402: res(3,node2) = res(3,node2) - res3
1403: res(4,node2) = res(4,node2) - res4
1404: flops = flops + 4.0
1405: endif
1406: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK)&& !defined(HAVE_EDGE_COLORING)
1407: endif
1408: #endif
1409: !
1410: enddo
1411: #if defined(_OPENMP)
1412: #if defined(HAVE_EDGE_COLORING)
1413: !$omp enddo
1414: !$omp end parallel
1415: first_edge = first_edge + ncount(i)
1416: enddo
1417: #elif defined(HAVE_REDUNDANT_WORK)
1418: !$omp enddo
1419: !$omp critical (residual_update)
1420: do n = 1,nnodes
1421: resvec(1,n)=resvec(1,n)+res(1,n)
1422: resvec(2,n)=resvec(2,n)+res(2,n)
1423: resvec(3,n)=resvec(3,n)+res(3,n)
1424: resvec(4,n)=resvec(4,n)+res(4,n)
1425: enddo
1426: flops = flops + 4.0*nnodes
1427: !$omp end critical (residual_update)
1428: !$omp end parallel
1429: #undef res
1430: #define res(a,b) resvec(a,b)
1431: #else
1432: !$omp end parallel
1433: #endif
1434: #else
1435: #endif
1436: call PetscLogFlops(flops,ierr)
1437: call PetscLogEventEnd(flux_event,ierr)
1438: ! c68 = 6./8.
1439: ! c18 = 1./8.
1440: c68 = 0.75d0
1441: c18 = 0.125d0
1442: !
1443: ! Close contour over the boundaries
1444: ! First do inviscid faces
1445: !
1446: do n = 1, nnfacet
1447: node1 = isnode(f2ntn(n,1))
1448: node2 = isnode(f2ntn(n,2))
1449: node3 = isnode(f2ntn(n,3))
1451: x1 = x(node1)
1452: y1 = y(node1)
1453: z1 = z(node1)
1454: p1 = qnod(1,node1)
1456: x2 = x(node2)
1457: y2 = y(node2)
1458: z2 = z(node2)
1459: p2 = qnod(1,node2)
1461: x3 = x(node3)
1462: y3 = y(node3)
1463: z3 = z(node3)
1464: p3 = qnod(1,node3)
1466: ax = x2 - x1
1467: ay = y2 - y1
1468: az = z2 - z1
1470: bx = x3 - x1
1471: by = y3 - y1
1472: bz = z3 - z1
1473: !
1474: ! Normal points away from grid interior.
1475: ! Magnitude is 1/3 area of surface triangle.
1476: !
1477: xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1478: ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1479: znorm =-0.5d0*(ax*by - ay*bx)/3.d0
1481: pa = c68*p1 + c18*(p2 + p3)
1482: pb = c68*p2 + c18*(p3 + p1)
1483: pc = c68*p3 + c18*(p1 + p2)
1484: !
1485: flops = flops + 35.0
1486: if (node1.le.nnodes) then
1487: res(2,node1) = res(2,node1) + xnorm*pa
1488: res(3,node1) = res(3,node1) + ynorm*pa
1489: res(4,node1) = res(4,node1) + znorm*pa
1490: flops = flops + 6.0
1491: endif
1493: if (node2.le.nnodes) then
1494: res(2,node2) = res(2,node2) + xnorm*pb
1495: res(3,node2) = res(3,node2) + ynorm*pb
1496: res(4,node2) = res(4,node2) + znorm*pb
1497: flops = flops + 6.0
1498: endif
1500: if (node3.le.nnodes) then
1501: res(2,node3) = res(2,node3) + xnorm*pc
1502: res(3,node3) = res(3,node3) + ynorm*pc
1503: res(4,node3) = res(4,node3) + znorm*pc
1504: flops = flops + 6.0
1505: endif
1507: enddo
1508: !
1509: ! Now viscous faces
1510: !
1511: do n = 1, nvfacet
1512: node1 = ivnode(f2ntv(n,1))
1513: node2 = ivnode(f2ntv(n,2))
1514: node3 = ivnode(f2ntv(n,3))
1516: x1 = x(node1)
1517: y1 = y(node1)
1518: z1 = z(node1)
1519: p1 = qnod(1,node1)
1521: x2 = x(node2)
1522: y2 = y(node2)
1523: z2 = z(node2)
1524: p2 = qnod(1,node2)
1526: x3 = x(node3)
1527: y3 = y(node3)
1528: z3 = z(node3)
1529: p3 = qnod(1,node3)
1531: ax = x2 - x1
1532: ay = y2 - y1
1533: az = z2 - z1
1535: bx = x3 - x1
1536: by = y3 - y1
1537: bz = z3 - z1
1538: !
1539: ! norm point away from grid interior.
1540: ! norm magnitude is 1/3 area of surface triangle.
1541: !
1542: xnorm =-0.5d0*(ay*bz - az*by)/3.d0
1543: ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
1544: znorm =-0.5d0*(ax*by - ay*bx)/3.d0
1546: pa = c68*p1 + c18*(p2 + p3)
1547: pb = c68*p2 + c18*(p3 + p1)
1548: pc = c68*p3 + c18*(p1 + p2)
1549: !
1550: flops = flops + 35.0
1551: !
1552: if (node1.le.nnodes) then
1553: res(2,node1) = res(2,node1) + xnorm*pa
1554: res(3,node1) = res(3,node1) + ynorm*pa
1555: res(4,node1) = res(4,node1) + znorm*pa
1556: flops = flops + 6.0
1558: endif
1560: if (node2.le.nnodes) then
1561: res(2,node2) = res(2,node2) + xnorm*pb
1562: res(3,node2) = res(3,node2) + ynorm*pb
1563: res(4,node2) = res(4,node2) + znorm*pb
1564: flops = flops + 6.0
1565: endif
1567: if (node3.le.nnodes) then
1568: res(2,node3) = res(2,node3) + xnorm*pc
1569: res(3,node3) = res(3,node3) + ynorm*pc
1570: res(4,node3) = res(4,node3) + znorm*pc
1571: flops = flops + 6.0
1572: endif
1574: enddo
1575: !
1576: ! The next section of code is for when you dont care about
1577: ! preserving linear data on boundary. Also, doing this
1578: ! matches the left hand side when not doing Newton-Krylov
1579: ! Usually just go around unless you are just experimenting
1580: !
1581: goto 1025
1582: !
1583: ! Loop over the boundaries
1584: ! First do inviscid nodes
1585: !
1586: ! do 1020 i = 1,nsnode
1587: ! inode = isnode(i)
1588: !
1589: ! xnorm = sxn(i)
1590: ! ynorm = syn(i)
1591: ! znorm = szn(i)
1592: !
1593: ! p = qnod(1,inode)
1594: !
1595: ! if (inode .le. nnodes) then
1596: ! res(2,inode) = res(2,inode) + xnorm*p
1597: ! res(3,inode) = res(3,inode) + ynorm*p
1598: ! res(4,inode) = res(4,inode) + znorm*p
1599: ! endif
1600: !
1601: !1020 continue
1602: !
1603: ! Now viscous nodes
1604: !
1605: ! do 1030 i = 1,nvnode
1606: ! inode = ivnode(i)
1607: !
1608: ! xnorm = vxn(i)
1609: ! ynorm = vyn(i)
1610: ! znorm = vzn(i)
1611: !
1612: ! p = qnod(1,inode)
1613: !
1614: ! if (inode .le. nnodes) then
1615: ! res(2,inode) = res(2,inode) + xnorm*p
1616: ! res(3,inode) = res(3,inode) + ynorm*p
1617: ! res(4,inode) = res(4,inode) + znorm*p
1618: ! endif
1619: !
1620: !1030 continue
1621: !
1622: 1025 continue
1623: !
1624: ! Now do far-field
1625: !
1626: do i = 1,nfnode
1627: inode = ifnode(i)
1628: !
1629: !/*
1630: ! Get normal and "other" 2 vectors. Remember that fxn,fyn and fzn
1631: ! has the magnitude of the face contained in it.
1632: !*/
1633: xnorm = fxn(i)
1634: ynorm = fyn(i)
1635: znorm = fzn(i)
1636: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
1637: xnorm = xnorm/area
1638: ynorm = ynorm/area
1639: znorm = znorm/area
1640: !
1641: ! Now lets get our other 2 vectors
1642: ! For first vector, use {1,0,0} and subtract off the component
1643: ! in the direction of the face normal. If the inner product of
1644: ! {1,0,0} is close to unity, use {0,1,0}
1645: !
1646: dot = xnorm
1647: if (abs(dot).lt.0.95d0)then
1648: X1 = 1.d0 - dot*xnorm
1649: Y1 = - dot*ynorm
1650: Z1 = - dot*znorm
1651: else
1652: dot = ynorm
1653: X1 = - dot*xnorm
1654: Y1 = 1.d0 - dot*ynorm
1655: Z1 = - dot*znorm
1656: end if
1657: !
1658: ! Normalize the first vector (V1)
1659: !
1660: size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
1661: X1 = X1/size
1662: Y1 = Y1/size
1663: Z1 = Z1/size
1664: !
1665: ! Take cross-product of normal with V1 to get V2
1666: !
1667: X2 = ynorm*Z1 - znorm*Y1
1668: Y2 = znorm*X1 - xnorm*Z1
1669: Z2 = xnorm*Y1 - ynorm*X1
1671: !
1672: ! Calculate elements of T and T(inverse) evaluated at freestream
1673: !
1674: ubar0 = xnorm*u0 + ynorm*v0 + znorm*w0
1675: c20 = ubar0*ubar0 + beta
1676: c0 = sqrt(c20)
1677: phi1 = xnorm*beta + u0*ubar0
1678: phi2 = ynorm*beta + v0*ubar0
1679: phi3 = znorm*beta + w0*ubar0
1680: phi4 = Y2*phi3 - Z2*phi2
1681: phi5 = Z2*phi1 - X2*phi3
1682: phi6 = X2*phi2 - Y2*phi1
1683: phi7 = Z1*phi2 - Y1*phi3
1684: phi8 = X1*phi3 - Z1*phi1
1685: phi9 = Y1*phi1 - X1*phi2
1687: t11 = 0.0d0
1688: t21 = X1
1689: t31 = Y1
1690: t41 = Z1
1692: t12 = 0.0d0
1693: t22 = X2
1694: t32 = Y2
1695: t42 = Z2
1697: t13 = c0*beta
1698: t23 = xnorm*beta + u0*(ubar0 + c0)
1699: t33 = ynorm*beta + v0*(ubar0 + c0)
1700: t43 = znorm*beta + w0*(ubar0 + c0)
1702: t14 = -c0*beta
1703: t24 = xnorm*beta + u0*(ubar0 - c0)
1704: t34 = ynorm*beta + v0*(ubar0 - c0)
1705: t44 = znorm*beta + w0*(ubar0 - c0)
1707: ti11 = -(u0*phi4 + v0*phi5 + w0*phi6)/beta
1708: ti21 = -(u0*phi7 + v0*phi8 + w0*phi9)/beta
1709: ti31 = .5d0*(c0 - ubar0)/beta
1710: ti41 = -.5d0*(c0 + ubar0)/beta
1711: !
1712: ! Now, get the variables on the "inside"
1713: !
1714: pi = qnod(1,inode)
1715: ui = qnod(2,inode)
1716: vi = qnod(3,inode)
1717: wi = qnod(4,inode)
1718: unorm = xnorm*ui + ynorm*vi + znorm*wi
1719: !
1720: ! If ubar is negative, take the reference condition from outside
1721: !
1722: !
1723: if (unorm.gt.0.0d0)then
1724: pr = pi
1725: ur = ui
1726: vr = vi
1727: wr = wi
1728: else
1729: pr = p0
1730: ur = u0
1731: vr = v0
1732: wr = w0
1733: end if
1734: !
1735: ! Set rhs
1736: !
1737: rhs1 = (ti11*pr + phi4*ur + phi5*vr + phi6*wr)/c20
1738: rhs2 = (ti21*pr + phi7*ur + phi8*vr + phi9*wr)/c20
1739: rhs3 = .5d0*(2.d0*ti31*pi + xnorm*ui + ynorm*vi + znorm*wi)/c20
1740: rhs4 = .5d0*(2.d0*ti41*p0 + xnorm*u0 + ynorm*v0 + znorm*w0)/c20
1741: !
1742: ! Now do matrix multiplication to get values on boundary
1743: !
1744: pb = t13*rhs3 + t14*rhs4
1745: ub = t21*rhs1 + t22*rhs2 + t23*rhs3 + t24*rhs4
1746: vb = t31*rhs1 + t32*rhs2 + t33*rhs3 + t34*rhs4
1747: wb = t41*rhs1 + t42*rhs2 + t43*rhs3 + t44*rhs4
1749: ubar = xnorm*ub + ynorm*vb + znorm*wb
1750: flops = flops + 180.0
1751: !
1752: if (inode.le.nnodes) then
1753: res(1,inode) = res(1,inode)+area*beta*ubar
1754: res(2,inode) = res(2,inode)+area*(ub*ubar + xnorm*pb)
1755: res(3,inode) = res(3,inode)+area*(vb*ubar + ynorm*pb)
1756: res(4,inode) = res(4,inode)+area*(wb*ubar + znorm*pb)
1757: flops = flops + 18.0
1758: endif
1759: !
1760: enddo
1761: call PetscLogFlops(flops,ierr)
1762: !
1763: ! End of subroutine FLUX
1764: !
1765: return
1766: end
1768: !---------------------------------------------------------------
1769: ! The following subroutines are from node3t.f in the original
1770: ! code - D. K. Kaushik (1/17/97)
1771: !---------------------------------------------------------------
1772: !
1773: !=============================== SUMGS ===============================72
1774: !
1775: ! Gets the weights for calculating gradients using least squares
1776: !
1777: !=====================================================================72
1778: subroutine SUMGS(nnodes,nedge,evec,xyz,rxy,irank,nvertices)
1780: implicit none
1781: #include <petsc/finclude/petscsys.h>
1782: integer nnodes,nedge,irank,nvertices
1783: #if defined(INTERLACING)
1784: PetscScalar xyz(3,nvertices)
1785: PetscScalar rxy(7,nnodes)
1786: integer evec(2,nedge)
1787: #define x(i) xyz(1,i)
1788: #define y(i) xyz(2,i)
1789: #define z(i) xyz(3,i)
1790: #define r11(i) rxy(1,i)
1791: #define r12(i) rxy(2,i)
1792: #define r13(i) rxy(3,i)
1793: #define r22(i) rxy(4,i)
1794: #define r23(i) rxy(5,i)
1795: #define r33(i) rxy(6,i)
1796: #define r44(i) rxy(7,i)
1797: #undef eptr
1798: #define eptr(j,i) evec(i,j)
1799: #else
1800: PetscScalar xyz(nvertices,3)
1801: PetscScalar rxy(nnodes,7)
1802: integer evec(nedge,2)
1803: #define x(i) xyz(i,1)
1804: #define y(i) xyz(i,2)
1805: #define z(i) xyz(i,3)
1806: #define r11(i) rxy(i,1)
1807: #define r12(i) rxy(i,2)
1808: #define r13(i) rxy(i,3)
1809: #define r22(i) rxy(i,4)
1810: #define r23(i) rxy(i,5)
1811: #define r33(i) rxy(i,6)
1812: #define r44(i) rxy(i,7)
1813: #define eptr(i,j) evec(i,j)
1814: #endif
1815: integer i,n,node1,node2,ierr
1816: PetscScalar x1,y1,z1,x2,y2,z2,dx,dy,dz
1817: PetscScalar dist,weight,w2,w11,w22,w33
1818: PetscScalar r12r11,r13r11,r23r22,rmess
1820: !
1821: ! Initialize all the rij to 0.0
1822: !
1823: do 1000 i = 1,nnodes
1824: r11(i) = 0.0
1825: r12(i) = 0.0
1826: r13(i) = 0.0
1827: r22(i) = 0.0
1828: r23(i) = 0.0
1829: r33(i) = 0.0
1830: r44(i) = 0.0
1831: 1000 continue
1832: !
1833: ! Now loop over the edges and accumulate the r
1834: !
1835: do 1020 n = 1, nedge
1836: !
1837: node1 = eptr(n,1)
1838: node2 = eptr(n,2)
1839: !
1840: x1 = x(node1)
1841: y1 = y(node1)
1842: z1 = z(node1)
1843: x2 = x(node2)
1844: y2 = y(node2)
1845: z2 = z(node2)
1846: !
1847: dx = x2 - x1
1848: dy = y2 - y1
1849: dz = z2 - z1
1850: dist = sqrt(dx*dx + dy*dy + dz*dz)
1851: ! weight = 1.0d0/dist
1852: weight = 1.0d0
1853: w2 = weight*weight
1854: !
1855: if (node1 .le. nnodes) then
1856: r11(node1) = r11(node1) + (x2 - x1)*(x2 - x1)*w2
1857: r12(node1) = r12(node1) + (x2 - x1)*(y2 - y1)*w2
1858: r13(node1) = r13(node1) + (x2 - x1)*(z2 - z1)*w2
1859: endif
1860: !
1861: if (node2 .le. nnodes) then
1862: r11(node2) = r11(node2) + (x1 - x2)*(x1 - x2)*w2
1863: r12(node2) = r12(node2) + (x1 - x2)*(y1 - y2)*w2
1864: r13(node2) = r13(node2) + (x1 - x2)*(z1 - z2)*w2
1865: endif
1866: 1020 continue
1867: !
1868: !/*
1869: ! Now calculate ||x(:)|| = r11 by taking the square root
1870: ! Also divide r12 and r13 by ||x(:)||
1871: !*/
1872: do 1030 i = 1,nnodes
1873: r11(i) = sqrt(r11(i))
1874: r12(i) = r12(i)/r11(i)
1875: r13(i) = r13(i)/r11(i)
1876: 1030 continue
1877: !
1878: !/*
1879: ! Now calculate r22 and r23
1880: !*/
1881: do 1050 n = 1, nedge
1882: !
1883: node1 = eptr(n,1)
1884: node2 = eptr(n,2)
1885: !
1886: x1 = x(node1)
1887: y1 = y(node1)
1888: z1 = z(node1)
1889: x2 = x(node2)
1890: y2 = y(node2)
1891: z2 = z(node2)
1892: !
1893: dx = x2 - x1
1894: dy = y2 - y1
1895: dz = z2 - z1
1896: dist = sqrt(dx*dx + dy*dy + dz*dz)
1897: ! weight = 1.0/dist
1898: weight = 1.0d0
1899: dx = weight*dx
1900: dy = weight*dy
1901: dz = weight*dz
1902: w2 = weight*weight
1903: !
1904: if (node1 .le. nnodes) then
1905: r22(node1) = r22(node1) + &
1906: & (dy-dx*r12(node1)/r11(node1))**2
1907: r23(node1) = r23(node1) + dz* &
1908: & (dy-dx*r12(node1)/r11(node1))
1909: endif
1910: !
1911: if (node2 .le. nnodes) then
1912: r22(node2) = r22(node2) + &
1913: & (-dy + dx*r12(node2)/r11(node2))**2
1914: r23(node2) = r23(node2) - dz* &
1915: & (-dy+dx*r12(node2)/r11(node2))
1916: endif
1917: 1050 continue
1918: !
1919: !/*
1920: ! Now finish getting r22 and r23
1921: !*/
1922: do 1060 i = 1,nnodes
1923: r22(i) = sqrt(r22(i))
1924: r23(i) = r23(i)/r22(i)
1925: 1060 continue
1926: !
1927: !/*
1928: ! Now all we have to do is get r33
1929: !*/
1930: do 1080 n = 1, nedge
1931: !
1932: node1 = eptr(n,1)
1933: node2 = eptr(n,2)
1934: !
1935: x1 = x(node1)
1936: y1 = y(node1)
1937: z1 = z(node1)
1938: x2 = x(node2)
1939: y2 = y(node2)
1940: z2 = z(node2)
1941: !
1942: dx = x2 - x1
1943: dy = y2 - y1
1944: dz = z2 - z1
1945: dist = sqrt(dx*dx + dy*dy + dz*dz)
1946: ! weight = 1.0/dist
1947: weight = 1.0d0
1948: dx = weight*dx
1949: dy = weight*dy
1950: dz = weight*dz
1951: w2 = weight*weight
1952: !
1953: if (node1 .le. nnodes) then
1954: r33(node1) = r33(node1) + &
1955: & (dz-dx*r13(node1)/r11(node1)- &
1956: & r23(node1)/r22(node1)* &
1957: & (dy - dx*r12(node1)/ &
1958: & r11(node1)))**2
1959: endif
1960: !
1961: if (node2 .le. nnodes) then
1962: r33(node2) = r33(node2) + &
1963: & (-dz+dx*r13(node2)/r11(node2)- &
1964: & r23(node2)/r22(node2)* &
1965: & (-dy + dx*r12(node2)/ &
1966: & r11(node2)))**2
1967: endif
1968: !
1969: 1080 continue
1970: !
1971: !/*
1972: ! Now just get the magnitude of r33
1973: !*/
1974: do 1090 i = 1,nnodes
1975: r33(i) = sqrt(r33(i))
1976: 1090 continue
1977: !
1978: !/*
1979: ! Added by Dinesh Kaushik (6/27/97)
1980: ! The following addition changes the meaning of r11 .. r33. r44
1981: ! is the new parameter introduced by me. The new definitions
1982: ! are taken from LSTGS (where these parameters
1983: ! are used finally).
1984: ! r11->w11
1985: ! r22->w22
1986: ! r33->w33
1987: ! r12->r12r11
1988: ! r13->r13r11
1989: ! r23->r23r22
1990: ! r44->rmess */
1992: do i = 1,nnodes
1993: w11 = 1.d0/(r11(i)*r11(i))
1994: w22 = 1.d0/(r22(i)*r22(i))
1995: w33 = 1.d0/(r33(i)*r33(i))
1996: r12r11 = r12(i)/r11(i)
1997: r13r11 = r13(i)/r11(i)
1998: r23r22 = r23(i)/r22(i)
1999: rmess=(r12(i)*r23(i)- &
2000: & r13(i)*r22(i))/ &
2001: & (r11(i)*r22(i)* &
2002: & r33(i)*r33(i))
2004: r11(i) = w11
2005: r22(i) = w22
2006: r33(i) = w33
2007: r12(i) = r12r11
2008: r13(i) = r13r11
2009: r23(i) = r23r22
2010: r44(i) = rmess
2011: enddo
2012: !
2013: ! Finished with SUMGS
2014: !
2015: return
2016: end
2018: !================================= LSTGS =============================72
2019: !
2020: ! Calculates the Gradients at the nodes using weighted least squares
2021: ! This subroutine solves using Gram-Schmidt
2022: !
2023: !=====================================================================72
2024: subroutine LSTGS(nnodes,nedge,evec, &
2025: & qvec,grad,xyz, &
2026: & rxy,irank,nvertices)
2027: !
2028: implicit none
2029: #include <petsc/finclude/petscsys.h>
2030: integer nnodes,nedge,irank,nvertices
2031: !
2032: #if defined(INTERLACING)
2033: PetscScalar qvec(4,nvertices)
2034: PetscScalar xyz(3,nvertices)
2035: PetscScalar rxy(7,nnodes)
2036: PetscScalar grad(3,4,nnodes)
2037: integer evec(2,nedge)
2038: #define qnod(i,j) qvec(i,j)
2039: #define gradx(i,j) grad(1,i,j)
2040: #define grady(i,j) grad(2,i,j)
2041: #define gradz(i,j) grad(3,i,j)
2042: #define x(i) xyz(1,i)
2043: #define y(i) xyz(2,i)
2044: #define z(i) xyz(3,i)
2045: #define r11(i) rxy(1,i)
2046: #define r12(i) rxy(2,i)
2047: #define r13(i) rxy(3,i)
2048: #define r22(i) rxy(4,i)
2049: #define r23(i) rxy(5,i)
2050: #define r33(i) rxy(6,i)
2051: #define r44(i) rxy(7,i)
2052: #define eptr(j,i) evec(i,j)
2053: #else
2054: PetscScalar qvec(nvertices,4)
2055: PetscScalar xyz(nvertices,3)
2056: PetscScalar rxy(nnodes,7)
2057: PetscScalar grad(nnodes,4,3)
2058: integer evec(nedge,2)
2059: #define qnod(i,j) qvec(j,i)
2060: #define gradx(i,j) grad(j,i,1)
2061: #define grady(i,j) grad(j,i,2)
2062: #define gradz(i,j) grad(j,i,3)
2063: #define x(i) xyz(i,1)
2064: #define y(i) xyz(i,2)
2065: #define z(i) xyz(i,3)
2066: #define r11(i) rxy(i,1)
2067: #define r12(i) rxy(i,2)
2068: #define r13(i) rxy(i,3)
2069: #define r22(i) rxy(i,4)
2070: #define r23(i) rxy(i,5)
2071: #define r33(i) rxy(i,6)
2072: #define r44(i) rxy(i,7)
2073: #define eptr(i,j) evec(i,j)
2074: #endif
2075: PetscScalar title(20),beta,alpha
2076: PetscScalar Re,dt,tot,res0,resc
2077: integer i,ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2078: PetscScalar cfl1,cfl2
2079: integer nsmoth,iflim,itran,nbtran,jupdate, &
2080: & nstage,ncyct,iramp,nitfo
2081: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
2082: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2083: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
2084: & nstage,ncyct,iramp,nitfo
2085: integer n,node1,node2,ierr
2086: PetscLogDouble flops
2087: PetscScalar dx1,dy1,dz1,dx2,dy2,dz2, &
2088: & dq1,dq2,dq3,dq4, &
2089: & weight,w11,r12r11,r13r11,w22,r23r22,w33,rmess, &
2090: & coef1,coef2,termx,termy,termz
2091: ! logging variables
2092: ! integer flag
2093: ! integer grad_event,node1_event,node2_event
2094: ! character * 16 grad_label, node1_label, node2_label
2095: ! data flag/-1/,grad_label/'GRAD '/
2096: ! data node1_label/'NODE1 '/
2097: ! data node2_label/'NODE2 '/
2098: ! save grad_event, grad_label, flag
2099: ! save node1_event,node2_event,node1_label, node2_label
2101: ! if (flag .eq. -1) then
2102: ! call PetscLogEventRegister(grad_event,grad_label,ierr)
2103: ! call PetscLogEventRegister(node1_event,node1_label,ierr)
2104: ! call PetscLogEventRegister(node2_event,node2_label,ierr)
2105: ! flag = 1
2106: ! endif
2107: ! call PetscLogEventBegin(grad_event,0,0,0,0,ierr)
2109: flops = 0.0
2110: ! For checking out the code input a linear distribution
2111: !
2112: ! write(6,700)nnodes,nedge
2113: ! 700 format(1h ,'nnodes=',i5,' nedge=',i5)
2114: ! do 1001 i = 1,nnodes
2115: ! write(6,800)i,x(i),y(i),z(i)
2116: ! 800 format(1h ,'i x y z=',i5,3(f10.5,1x))
2117: ! qnod(1,i) = 1.0*x(i) + 2.0*y(i) + 3.0*z(i)
2118: ! qnod(2,i) = 3.0*x(i) + 4.0*y(i) + 6.0*z(i)
2119: ! qnod(3,i) = 5.0*x(i) + 6.0*y(i) + 9.0*z(i)
2120: ! qnod(4,i) = 7.0*x(i) + 8.0*y(i) + 12.0*z(i)
2121: ! qnod(i,5) = 9.0*x(i) + 10.0*y(i) + 15.0*z(i)
2122: !1001 continue
2123: !
2124: ! Zero out the gradients
2125: !
2126: do 1000 i = 1,nnodes
2127: gradx(1,i) = 0.0
2128: grady(1,i) = 0.0
2129: gradz(1,i) = 0.0
2131: gradx(2,i) = 0.0
2132: grady(2,i) = 0.0
2133: gradz(2,i) = 0.0
2135: gradx(3,i) = 0.0
2136: grady(3,i) = 0.0
2137: gradz(3,i) = 0.0
2139: gradx(4,i) = 0.0
2140: grady(4,i) = 0.0
2141: gradz(4,i) = 0.0
2142: 1000 continue
2144: ! If second order, loop over all the faces accumulate sums
2145: !
2146: ! nitfo = 0
2147: ! if (ntt.gt.nitfo.or.ivisc.gt.0)then
2148: ! if (1 .gt. 0) then
2149: do 1020 n = 1, nedge
2150: node1 = eptr(n,1)
2151: node2 = eptr(n,2)
2152: ! if ((node1 .le. nnodes).or.(node2 .le. nnodes)) then
2153: dx1 = x(node2) - x(node1)
2154: dy1 = y(node2) - y(node1)
2155: dz1 = z(node2) - z(node1)
2156: !
2157: ! dist = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1)
2158: ! weight = 1.0/dist
2159: weight = 1.0d0
2160: ! w2 = weight*weight
2161: !
2162: dx1 = weight*dx1
2163: dy1 = weight*dy1
2164: dz1 = weight*dz1
2166: flops = flops + 6.0
2167: !
2168: ! call PetscLogEventBegin(node1_event,0,0,0,0,ierr)
2169: if (node1 .le. nnodes) then
2170: dq1 = weight*(qnod(1,node2) - qnod(1,node1))
2171: dq2 = weight*(qnod(2,node2) - qnod(2,node1))
2172: dq3 = weight*(qnod(3,node2) - qnod(3,node1))
2173: dq4 = weight*(qnod(4,node2) - qnod(4,node1))
2174: !
2175: ! w11 = 1./(r11(node1)*r11(node1))
2176: ! w22 = 1./(r22(node1)*r22(node1))
2177: ! w33 = 1./(r33(node1)*r33(node1))
2178: ! r12r11 = r12(node1)/r11(node1)
2179: ! r13r11 = r13(node1)/r11(node1)
2180: ! r23r22 = r23(node1)/r22(node1)
2181: ! rmess = (r12(node1)*r23(node1) - r13(node1)*r22(node1))/
2182: ! 1 (r11(node1)*r22(node1)*r33(node1)*r33(node1))
2183: !
2184: w11 = r11(node1)
2185: r12r11 = r12(node1)
2186: r13r11 = r13(node1)
2187: w22 = r22(node1)
2188: r23r22 = r23(node1)
2189: w33 = r33(node1)
2190: rmess = r44(node1)
2191: !
2192: coef1 = dy1 - dx1*r12r11
2193: coef2 = dz1 - dx1*r13r11 - r23r22*coef1
2194: termx = dx1*w11 - w22*r12r11*coef1 + rmess*coef2
2195: termy = w22*coef1 - r23r22*w33*coef2
2196: termz = w33*coef2
2197: !
2198: gradx(1,node1) = gradx(1,node1) + termx*dq1
2199: grady(1,node1) = grady(1,node1) + termy*dq1
2200: gradz(1,node1) = gradz(1,node1) + termz*dq1
2201: !
2202: gradx(2,node1) = gradx(2,node1) + termx*dq2
2203: grady(2,node1) = grady(2,node1) + termy*dq2
2204: gradz(2,node1) = gradz(2,node1) + termz*dq2
2205: !
2206: gradx(3,node1) = gradx(3,node1) + termx*dq3
2207: grady(3,node1) = grady(3,node1) + termy*dq3
2208: gradz(3,node1) = gradz(3,node1) + termz*dq3
2209: !
2210: gradx(4,node1) = gradx(4,node1) + termx*dq4
2211: grady(4,node1) = grady(4,node1) + termy*dq4
2212: gradz(4,node1) = gradz(4,node1) + termz*dq4
2213: !
2214: flops = flops + 49.0
2215: endif
2216: ! call PetscLogEventEnd(node1_event,0,0,0,0,ierr)
2217: !
2218: ! Now do the other node
2219: !
2220: ! call PetscLogEventBegin(node2_event,0,0,0,0,ierr)
2221: if (node2 .le. nnodes) then
2222: dx2 = -dx1
2223: dy2 = -dy1
2224: dz2 = -dz1
2225: !
2226: dq1 = weight*(qnod(1,node1) - qnod(1,node2))
2227: dq2 = weight*(qnod(2,node1) - qnod(2,node2))
2228: dq3 = weight*(qnod(3,node1) - qnod(3,node2))
2229: dq4 = weight*(qnod(4,node1) - qnod(4,node2))
2230: !
2231: ! w11 = 1./(r11(node2)*r11(node2))
2232: ! w22 = 1./(r22(node2)*r22(node2))
2233: ! w33 = 1./(r33(node2)*r33(node2))
2234: ! r12r11 = r12(node2)/r11(node2)
2235: ! r13r11 = r13(node2)/r11(node2)
2236: ! r23r22 = r23(node2)/r22(node2)
2237: ! rmess = (r12(node2)*r23(node2) - r13(node2)*r22(node2))/
2238: ! & (r11(node2)*r22(node2)*r33(node2)*r33(node2))
2240: w11 = r11(node2)
2241: r12r11 = r12(node2)
2242: r13r11 = r13(node2)
2243: w22 = r22(node2)
2244: r23r22 = r23(node2)
2245: w33 = r33(node2)
2246: rmess = r44(node2)
2247: !
2248: coef1 = dy2 - dx2*r12r11
2249: coef2 = dz2 - dx2*r13r11 - r23r22*coef1
2250: termx = dx2*w11 - w22*r12r11*coef1 + rmess*coef2
2251: termy = w22*coef1 - r23r22*w33*coef2
2252: termz = w33*coef2
2253: !
2254: gradx(1,node2) = gradx(1,node2) + termx*dq1
2255: grady(1,node2) = grady(1,node2) + termy*dq1
2256: gradz(1,node2) = gradz(1,node2) + termz*dq1
2257: !
2258: gradx(2,node2) = gradx(2,node2) + termx*dq2
2259: grady(2,node2) = grady(2,node2) + termy*dq2
2260: gradz(2,node2) = gradz(2,node2) + termz*dq2
2261: !
2262: gradx(3,node2) = gradx(3,node2) + termx*dq3
2263: grady(3,node2) = grady(3,node2) + termy*dq3
2264: gradz(3,node2) = gradz(3,node2) + termz*dq3
2265: !
2266: gradx(4,node2) = gradx(4,node2) + termx*dq4
2267: grady(4,node2) = grady(4,node2) + termy*dq4
2268: gradz(4,node2) = gradz(4,node2) + termz*dq4
2269: !
2270: flops = flops + 52.0
2271: endif
2272: ! call PetscLogEventEnd(node2_event,0,0,0,0,ierr)
2273: ! endif
2274: !
2275: 1020 continue
2276: ! end if
2277: call PetscLogFlops(flops,ierr)
2278: ! call PetscLogEventEnd(grad_event,0,0,0,0,ierr)
2279: !
2280: ! End of LSTGS
2281: !
2282: return
2283: end
2285: !=================================== GETRES ==========================72
2286: !
2287: ! Calculates the residual
2288: ! Last Modified - D. K. Kaushik 1/23/97
2289: ! I have eliminated the input variables which were not needed -
2290: ! dq, A, B, iupdate
2291: !
2292: !=====================================================================72
2293: #if defined(_OPENMP)
2294: #if defined(HAVE_EDGE_COLORING)
2295: subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface, &
2296: & nsnode,nvnode,nfnode,isface,ivface,ifface, &
2297: & ileast,isnode,ivnode,ifnode, &
2298: & nnfacet,f2ntn,nnbound, &
2299: & nvfacet,f2ntv,nvbound, &
2300: & nffacet,f2ntf,nfbound, &
2301: & evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
2302: & xyzn,qvec,cdt,xyz,area,grad,resvec,turbre, &
2303: & slen,c2n,c2e,us,vs,as,phi,amut,ires, &
2304: & max_threads, &
2305: & ncolor,ncount, &
2306: & LocalTS,irank, nvertices)
2307: #elif defined(HAVE_REDUNDANT_WORK)
2308: subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface, &
2309: & nsnode,nvnode,nfnode,isface,ivface,ifface, &
2310: & ileast,isnode,ivnode,ifnode, &
2311: & nnfacet,f2ntn,nnbound, &
2312: & nvfacet,f2ntv,nvbound, &
2313: & nffacet,f2ntf,nfbound, &
2314: & evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
2315: & xyzn,qvec,cdt,xyz,area,grad,resvec,turbre, &
2316: & slen,c2n,c2e,us,vs,as,phi,amut,ires, &
2317: & max_threads, &
2318: & resd, &
2319: & LocalTS,irank, nvertices)
2320: #else
2321: subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface, &
2322: & nsnode,nvnode,nfnode,isface,ivface,ifface, &
2323: & ileast,isnode,ivnode,ifnode, &
2324: & nnfacet,f2ntn,nnbound, &
2325: & nvfacet,f2ntv,nvbound, &
2326: & nffacet,f2ntf,nfbound, &
2327: & evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
2328: & xyzn,qvec,cdt,xyz,area,grad,resvec,turbre, &
2329: & slen,c2n,c2e,us,vs,as,phi,amut,ires, &
2330: & max_threads, &
2331: & nedgeAllThr,part_thr,nedge_thr, &
2332: & edge_thr,xyzn_thr, &
2333: & LocalTS,irank, nvertices)
2334: #endif
2335: #else
2336: subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface, &
2337: & nsnode,nvnode,nfnode,isface,ivface,ifface, &
2338: & ileast,isnode,ivnode,ifnode, &
2339: & nnfacet,f2ntn,nnbound, &
2340: & nvfacet,f2ntv,nvbound, &
2341: & nffacet,f2ntf,nfbound, &
2342: & evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
2343: & xyzn,qvec,cdt,xyz,area,grad,resvec,turbre, &
2344: & slen,c2n,c2e,us,vs,as,phi,amut,ires, &
2345: & LocalTS,irank, nvertices)
2346: #endif
2347: !
2348: implicit none
2349: #include <petsc/finclude/petscsys.h>
2350: !
2351: integer nnodes, ncell, nedge, &
2352: & nbface,ileast,ires, &
2353: & nsface, nvface, nfface, &
2354: & nsnode, nvnode, nfnode, &
2355: & nnfacet,nvfacet,nffacet, &
2356: & nnbound,nvbound,nfbound, &
2357: & LocalTS,irank,nvertices
2358: integer evec(nedge,2)
2359: integer isface(1),ivface(1),ifface(1)
2360: integer isnode(1),ivnode(1),ifnode(1)
2361: integer c2n(ncell,4),c2e(ncell,6)
2362: integer f2ntn(nnfacet,4)
2363: integer f2ntv(nvfacet,4)
2364: integer f2ntf(nffacet,4)
2365: !
2366: PetscScalar us(nbface,3,4),vs(nbface,3,4)
2367: PetscScalar as(nbface,3,4)
2368: PetscScalar sxn(1),syn(1),szn(1)
2369: PetscScalar vxn(1),vyn(1),vzn(1)
2370: PetscScalar fxn(1),fyn(1),fzn(1)
2371: PetscScalar xyz(nvertices,3),area(nvertices)
2372: PetscScalar xyzn(nedge,4)
2373: PetscScalar turbre(1),slen(1)
2374: PetscScalar qvec(nvertices,4)
2375: PetscScalar phi(nvertices,4)
2376: PetscScalar cdt(nvertices)
2377: PetscScalar grad(3,4,nvertices)
2378: ! real dq(nnodes,4)
2379: ! real r11(nvertices),r12(nvertices),r13(nvertices)
2380: ! real r22(nvertices),r23(nvertices),r33(nvertices)
2381: PetscScalar amut(nnodes)
2382: !
2383: integer i
2384: PetscScalar title(20),beta,alpha
2385: PetscScalar Re,dt,tot,res0,resc
2386: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2387: PetscScalar cfl1,cfl2
2388: integer nsmoth,iflim,itran,nbtran,jupdate, &
2389: & nstage,ncyct,iramp,nitfo
2390: PetscScalar gtol
2391: integer icycle,nsrch,ilu0,ifcn
2392: PetscScalar rms(1001),clw(1001),cdw(1001), &
2393: & cmw(1001),xres(1001)
2394: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
2395: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
2396: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
2397: & nstage,ncyct,iramp,nitfo
2398: common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
2399: common/history/rms,clw,cdw,cmw,xres
2400: ! logging variables
2401: ! integer flux_event, flag, ierr
2402: ! integer delta2_event
2403: ! character * 16 flux_label
2404: ! character * 16 delta2_label
2405: ! data flag/-1/,flux_label/'FluxEvaluation '/
2406: ! data delta2_label/'DELTA2 '/
2407: ! save flux_event, flag, flux_label
2408: ! save delta2_event,delta2_label
2409: #if defined(INTERLACING)
2410: PetscScalar resvec(4,nnodes)
2411: #undef res
2412: #define res(i,j) resvec(i,j)
2413: #else
2414: PetscScalar resvec(nnodes,4)
2415: #define res(i,j) resvec(j,i)
2416: #endif
2417: #if defined(_OPENMP)
2418: integer max_threads
2419: #if defined(HAVE_EDGE_COLORING)
2420: integer ncolor,ncount(1)
2421: #elif defined(HAVE_REDUNDANT_WORK)
2422: PetscScalar resd(4,nnodes)
2423: #else
2424: integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads), &
2425: & edge_thr(2,nedgeAllThr)
2426: PetscScalar xyzn_thr(4,nedgeAllThr)
2427: #endif
2428: #endif
2429: !
2431: ! if (flag .eq. -1) then
2432: ! call PetscLogEventRegister(delta2_event,delta2_label,ierr)
2433: ! call PetscLogEventRegister(flux_event,flux_label,ierr)
2434: ! flag = 1
2435: ! endif
2436: !
2437: ! Calculate the time step
2438: !
2439: if (ires.eq.1) goto 888
2440: ! call PetscLogEventBegin(delta2_event,0,0,0,0,ierr)
2441: call DELTAT2(nnodes,nedge,qvec,cdt, &
2442: & xyz,area,xyzn,evec, &
2443: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn, &
2444: & nsnode,nvnode,nfnode,isnode,ivnode,ifnode, &
2445: & LocalTS,irank,nvertices)
2446: ! call PetscLogEventEnd(delta2_event,0,0,0,0,ierr)
2447: 888 continue
2448: !
2449: !/*
2450: ! Calculate the gradients
2451: ! ----Kyle seems to recommend only LSTGS for gradients,
2452: ! so I have commented the GETGRAD call - DKK (1/17/97)
2453: !
2454: ! if (ileast.eq.0) then
2455: ! call GETGRAD(nnodes,ncell,nedge,nsface,nvface,nfface,
2456: ! & isface,ivface,ifface,eptr,ncolor,ncount,
2457: ! & qnod,gradx,grady,x,y,
2458: ! & area,wx,wy,xn,yn,rl)
2459: !
2460: ! else if (ileast.eq.4) then
2461: ! if (ileast.eq.4) then
2462: ! call LSTGS(nnodes,nedge,eptr,
2463: ! & qnod,gradx,grady,gradz,xyz,
2464: ! & r11,r12,r13,r22,r23,r33,irank,nvertices)
2465: ! end if
2466: !
2467: ! zero out residuals (viscous residuals are zeroed in vfluxnew)
2468: !*/
2469: do 1002 i = 1,nnodes
2470: res(1,i)=0.0d0
2471: res(2,i)=0.0d0
2472: res(3,i)=0.0d0
2473: res(4,i)=0.0d0
2475: phi(i,1)=1.0d0
2476: phi(i,2)=1.0d0
2477: phi(i,3)=1.0d0
2478: phi(i,4)=1.0d0
2479: 1002 continue
2480: !
2481: !/*
2482: !
2483: ! If not doing Newton-Krylov and iflim=1 call the Flux Limiter
2484: !
2485: ! if (iflim.eq.1.and.ifcn.ne.1)then
2486: ! call TIMLIM(nnodes,nedge,qnod,res,dq,phi,ncolor,ncount,
2487: ! & gradx,grady,gradz,x,y,z,eptr)
2488: !
2489: ! If we used the limiter we need to zero out the residual again
2490: ! since we used it for scratch space
2491: !
2492: ! do 1003 i = 1,nnodes
2493: ! res(1,i)=0.0
2494: ! res(2,i)=0.0
2495: ! res(3,i)=0.0
2496: ! res(4,i)=0.0
2497: !1003 continue
2498: !
2499: ! end if
2500: !*/
2501: ! Split the fluxes and perform the flux balance
2502: !
2503: ! call PetscLogEventBegin(flux_event,0,0,0,0,ierr)
2505: #if defined(_OPENMP)
2506: #if defined(HAVE_EDGE_COLORING)
2507: call FLUX(nnodes, ncell, nedge, &
2508: & nsface, nvface, nfface, isface, ivface, ifface, &
2509: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
2510: & nnfacet,f2ntn,nnbound, &
2511: & nvfacet,f2ntv,nvbound, &
2512: & nffacet,f2ntf,nfbound, &
2513: & grad,evec,qvec,xyz,resvec,xyzn, &
2514: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
2515: & max_threads, &
2516: & ncolor,ncount, &
2517: & irank,nvertices)
2518: #elif defined(HAVE_REDUNDANT_WORK)
2519: call FLUX(nnodes, ncell, nedge, &
2520: & nsface, nvface, nfface, isface, ivface, ifface, &
2521: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
2522: & nnfacet,f2ntn,nnbound, &
2523: & nvfacet,f2ntv,nvbound, &
2524: & nffacet,f2ntf,nfbound, &
2525: & grad,evec,qvec,xyz,resvec,xyzn, &
2526: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
2527: & max_threads, &
2528: & resd, &
2529: & irank,nvertices)
2530: #else
2531: call FLUX(nnodes, ncell, nedge, &
2532: & nsface, nvface, nfface, isface, ivface, ifface, &
2533: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
2534: & nnfacet,f2ntn,nnbound, &
2535: & nvfacet,f2ntv,nvbound, &
2536: & nffacet,f2ntf,nfbound, &
2537: & grad,evec,qvec,xyz,resvec,xyzn, &
2538: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
2539: & max_threads, &
2540: & nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr, &
2541: & irank,nvertices)
2542: #endif
2543: #else
2544: call FLUX(nnodes, ncell, nedge, &
2545: & nsface, nvface, nfface, isface, ivface, ifface, &
2546: & nsnode, nvnode, nfnode, isnode, ivnode, ifnode, &
2547: & nnfacet,f2ntn,nnbound, &
2548: & nvfacet,f2ntv,nvbound, &
2549: & nffacet,f2ntf,nfbound, &
2550: & grad,evec,qvec,xyz,resvec,xyzn, &
2551: & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi, &
2552: & irank,nvertices)
2553: #endif
2554: ! call PetscLogEventEnd(flux_event,0,0,0,0,ierr)
2555: !/*
2556: ! calculate viscous fluxes
2557: !
2558: ! if (ivisc.gt.0) then
2559: ! call VISRHS (nnodes,ncell,nedge,
2560: ! call EDGEVIS(nnodes,ncell,nedge,
2561: ! & nsnode,nvnode,nfnode,isnode,ivnode,ifnode,
2562: ! & nsface,nvface,nfface,isface,ivface,ifface,
2563: ! & nnfacet,f2ntn,nnbound,ncolorn,countn,
2564: ! & nvfacet,f2ntv,nvbound,ncolorv,countv,
2565: ! & nffacet,f2ntf,nfbound,ncolorf,countf,
2566: ! & nccolor,nccount,
2567: ! & eptr,c2n,c2e,
2568: ! & sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,
2569: ! & x,y,z,gradx,grady,gradz,
2570: ! & qnod,amut,res,phi)
2571: ! end if
2572: !*/
2573: ! End of subroutine GETRES
2574: !
2575: return
2576: end
2578: !---------------------------------------------------------------
2579: ! The following subroutine is from node4t.f in the original
2580: ! code - D. K. Kaushik (1/17/97)
2581: !---------------------------------------------------------------
2582: !
2583: !=============================================================================
2584: !
2585: ! Opens files for I/O
2586: !
2587: !=============================================================================
2588: SUBROUTINE OPENM(irank)
2589: !
2590: ! TAPE7 -- input: mach number, angle of attack etc..
2591: ! TAPE9 -- input: reads restart file
2592: ! TAPE10 -- output: residual history
2593: ! TAPE11 -- output: writes restart file
2594: ! TAPE12 -- output: writes residual and lift for plotting
2595: ! TAPE13 -- output: writes flowfield for contour plotting
2596: !
2597: implicit none
2598: integer irank
2599: ! OPEN(UNIT= 7,FILE='home/petsc/datafiles/fun3dgrid/ginput.faces', &
2600: ! &form='formatted',STATUS='OLD')
2602: ! OPEN(UNIT=9,FILE='framer.bin',
2603: ! &form='unformatted',STATUS='old')
2605: if (irank .eq. 0) OPEN(UNIT= 10,FILE='frame.out3', &
2606: &form='formatted',STATUS='unknown')
2608: ! OPEN(UNIT= 11,FILE='frame.bin',
2609: ! +form='unformatted',STATUS='unknown')
2611: ! OPEN(UNIT= 12,FILE='frame.plt',
2612: ! +form='formatted',STATUS='unknown')
2614: ! OPEN(UNIT= 13,FILE='frame.tec',
2615: ! +form='formatted',STATUS='unknown')
2617: ! OPEN(UNIT= 14,FILE='frame.fast.g',
2618: ! +form='unformatted',STATUS='unknown')
2620: ! OPEN(UNIT= 15,FILE='frame.fast.q',
2621: ! +form='unformatted',STATUS='unknown')
2623: return
2624: end
2625: !
2626: !
2627: !===================================================================
2628: !
2629: ! Get the IA, JA, and IAU arrays
2630: !
2631: !===================================================================
2632: subroutine GETIA(nnodes,nedge,evec,ia,ideg,irank)
2633: implicit none
2634: #include <petsc/finclude/petscsys.h>
2635: integer nnodes,nedge,irank
2636: integer ia(1),ideg(1)
2637: #if defined(INTERLACING)
2638: integer evec(2,nedge)
2639: #define eptr(j,i) evec(i,j)
2640: #else
2641: integer evec(nedge,2)
2642: #define eptr(i,j) evec(i,j)
2643: #endif
2644: integer i,node1,node2
2645: !
2646: ! First get the degree of each node using ideg as a dummy array
2647: !
2648: do 1000 i = 1,nnodes
2649: ideg(i) = 0
2650: 1000 continue
2651: !
2652: do 1010 i = 1,nedge
2653: node1 = eptr(i,1)
2654: node2 = eptr(i,2)
2655: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
2656: ideg(node1) = ideg(node1) + 1
2657: ideg(node2) = ideg(node2) + 1
2658: #else
2659: if (node1 .le. nnodes) ideg(node1) = ideg(node1) + 1
2660: if (node2 .le. nnodes) ideg(node2) = ideg(node2) + 1
2661: #endif
2662: 1010 continue
2663: !
2664: ! Now we can fill the ia array fairly easily
2665: !
2666: ia(1) = 1
2667: do 1020 i = 1,nnodes
2668: ia(i+1) = ia(i) + ideg(i) + 1
2669: ! write(9,100)i,ideg(i)
2670: ! 100 format(1h ,'deg(',i6,')=',i6)
2671: 1020 continue
2672: !
2673: return
2674: end
2675: !
2676: !===================================================================
2677: !
2678: ! Get the IA, JA, and IAU arrays
2679: !
2680: !===================================================================
2681: subroutine GETJA(nnodes,nedge,evec,ia,ja,iwork,irank)
2682: implicit none
2683: #include <petsc/finclude/petscsys.h>
2684: integer nnodes,nedge,irank
2685: integer ia(1),ja(1),iwork(1)
2686: #if defined(INTERLACING)
2687: integer evec(2,nedge)
2688: #define eptr(j,i) evec(i,j)
2689: #else
2690: integer evec(nedge,2)
2691: #define eptr(i,j) evec(i,j)
2692: #endif
2693: integer i,index,node1,node2,index1,index2, &
2694: & istart,iend
2695: ! open(unit=90,file='map.out',status='UNKNOWN')
2696: !
2697: ! Now we need to get the JA array
2698: ! First fill the diagonal places
2699: !
2700: do 1040 i = 1,nnodes
2701: index = ia(i)
2702: ja(index) = i
2703: iwork(i) = 1
2704: 1040 continue
2705: !
2706: do 1030 i = 1,nedge
2707: node1 = eptr(i,1)
2708: node2 = eptr(i,2)
2709: !
2710: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
2711: index1 = ia(node1) + iwork(node1)
2712: iwork(node1) = iwork(node1) + 1
2713: ja(index1) = node2
2714: index2 = ia(node2) + iwork(node2)
2715: iwork(node2) = iwork(node2) + 1
2716: ja(index2) = node1
2717: #else
2718: if (node1 .le. nnodes) then
2719: index1 = ia(node1) + iwork(node1)
2720: iwork(node1) = iwork(node1) + 1
2721: ja(index1) = node2
2722: endif
2723: if (node2 .le. nnodes) then
2724: index2 = ia(node2) + iwork(node2)
2725: iwork(node2) = iwork(node2) + 1
2726: ja(index2) = node1
2727: endif
2728: #endif
2729: 1030 continue
2730: !
2731: ! Now lets sort all our "bins" and get the correct one on the diagonal
2732: !
2733: do 1050 i = 1,nnodes
2734: istart = ia(i)
2735: iend = ia(i+1) - 1
2736: ! write(9,200)i,istart,iend
2737: ! 200 format(1h ,'Sorting ',i6,' istart iend = ',i6,1x,i6)
2738: call SORTER(istart,iend,ja,i)
2739: 1050 continue
2740: !
2741: ! Now get the "fhelp" array which will assist in assembling
2742: ! the flux Jacobians into the correct location in the alu array
2743: !
2744: ! write(90,*) 'fhelp array'
2745: ! do 1060 i = 1,nedge
2746: ! node1 = eptr(i,1)
2747: ! node2 = eptr(i,2)
2748: !
2749: ! First take care of node1
2750: !
2751: ! idiag = iau(node1)
2752: !
2753: ! If the offdiagonal term is ordered later in the ja array
2754: !
2755: ! if (node2.gt.node1)then
2756: ! jstart = idiag + 1
2757: ! jend = ia(node1+1) - 1
2758: ! else
2759: ! jstart = ia(node1)
2760: ! jend = idiag -1
2761: ! end if
2762: !
2763: ! do 1070 j = jstart,jend
2764: ! if (ja(j).eq.node2)fhelp(i,1) = j
2765: !1070 continue
2766: !
2767: !
2768: ! Now take care of node2
2769: !
2770: ! idiag = iau(node2)
2771: !
2772: ! If the offdiagonal term is ordered later in the ja array
2773: !
2774: ! if (node1.gt.node2)then
2775: ! jstart = idiag + 1
2776: ! jend = ia(node2+1) - 1
2777: ! else
2778: ! jstart = ia(node2)
2779: ! jend = idiag -1
2780: ! end if
2781: !
2782: ! do 1080 j = jstart,jend
2783: ! if (ja(j).eq.node1)fhelp(i,2) = j
2784: !1080 continue
2785: ! write(90,*) i,fhelp(i,1),fhelp(i,2)
2786: !1060 continue
2787: ! close(90)
2788: !
2789: return
2790: end
2791: !
2792: !
2793: !===================================================================
2794: !
2795: ! Sort each of our bins
2796: !
2797: !===================================================================
2798: subroutine SORTER(istart,iend,ja,inode)
2799: implicit none
2800: #include <petsc/finclude/petscsys.h>
2801: integer istart,iend,inode
2802: integer ja(1)
2803: !
2804: integer min,minsave,jsave,i,j
2805: do 1000 i = istart,iend
2806: min = ja(i)
2807: minsave = ja(i)
2808: jsave = i
2809: do 1010 j = i+1,iend
2810: if (ja(j).lt.min)then
2811: min = ja(j)
2812: jsave = j
2813: end if
2814: 1010 continue
2815: ja(i) = min
2816: ja(jsave) = minsave
2817: ! if (ja(i).eq.inode)iau(inode) = i
2818: 1000 continue
2819: !
2820: return
2821: end