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