Actual source code: userJac.F
1: !---------------------------------------------------------------
2: ! The following subroutines are from node6t.f in the original
3: ! code.
4: ! Last Modified - D. K. Kaushik (1/24/97)
5: !---------------------------------------------------------------
6: !
7: !=================================== FILLA ===================================
8: !
9: ! Fill the nonzero term of the A matrix
10: ! Modified - D. K. Kaushik (1/14/97)
11: ! cfl is passed as a parameter
12: !
13: !=============================================================================
14: subroutine FILLA(nnodes,nedge,evec, &
15: & nsface,isface,fxn,fyn,fzn,sxn,syn,szn, &
16: & nsnode,nvnode,nfnode,isnode,ivnode,ifnode, &
17: & qvec,A,cdt, &
18: & vol,xyzn,cfl,irank,nvertices)
19: !
20: implicit none
21: #include <petsc/finclude/petscsys.h>
22: #include <petsc/finclude/petscvec.h>
23: #include <petsc/finclude/petscmat.h>
24: #include <petsc/finclude/petscpc.h>
25: #include <petsc/finclude/petscsnes.h>
26: integer nnodes,nedge,nsface,nsnode,nvnode,nfnode,irank,nvertices
28: integer isface(1)
29: integer isnode(1),ivnode(1),ifnode(1)
30: PetscScalar cfl
31: #if defined(INTERLACING)
32: PetscScalar qvec(4,nvertices)
33: PetscScalar xyzn(4,nedge)
34: integer evec(2,nedge)
35: #define qnode(i,j) qvec(i,j)
36: #define dfp(a,b) val1(b,a)
37: #define dfm(a,b) val1(b+4,a)
38: #define xn(i) xyzn(1,i)
39: #define yn(i) xyzn(2,i)
40: #define zn(i) xyzn(3,i)
41: #define rl(i) xyzn(4,i)
42: #define eptr(j,i) evec(i,j)
43: #else
44: PetscScalar qvec(nvertices,4)
45: PetscScalar xyzn(nedge,4)
46: integer evec(nedge,2)
47: #define qnode(i,j) qvec(j,i)
48: #define dfp(a,b) val1(b,a)
49: #define dfm(a,b) val1(b+4,a)
50: #define xn(i) xyzn(i,1)
51: #define yn(i) xyzn(i,2)
52: #define zn(i) xyzn(i,3)
53: #define rl(i) xyzn(i,4)
54: #define eptr(i,j) evec(i,j)
55: #endif
56: PetscScalar vol(1)
57: PetscScalar sxn(1),syn(1),szn(1)
58: PetscScalar fxn(1),fyn(1),fzn(1)
59: PetscScalar cdt(1)
60: ! PetscScalar A(nnz,4,4)
61: Mat A
62: MatStructure flag
64: ! integer ia(1),ja(1),iau(1),fhelp(nedge,2)
65: PetscScalar title(20),beta,alpha
66: PetscScalar Re,dt,tot,res0,resc
67: integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
68: PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
69: PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
70: PetscScalar cfl1,cfl2
71: integer nsmoth,iflim,itran,nbtran,jupdate,nstage,ncyct,iramp,nitfo
72: common/info/title,beta,alpha,Re,dt,tot,res0,resc, &
73: & ntt,mseq,ivisc,irest,icyc,ihane,ntturb
74: common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
75: common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
76: common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate, &
77: & nstage,ncyct,iramp,nitfo
78: integer irow(16), icol(16)
79: integer i,j,k,n,ic,ierr,ir,node1,node2,inode
80: PetscLogDouble flops
81: PetscScalar val(32), val1(8,4)
82: PetscScalar xnorm,ynorm,znorm,rlen,dot,temp
83: PetscScalar X1,Y1,Z1,X2,Y2,Z2,size,area
84: PetscScalar pL,uL,vL,wL,ubarL,c2L,cL
85: PetscScalar pR,uR,vR,wR,ubarR,c2R,cR
86: PetscScalar pi,ui,vi,wi,unorm,c20,ubar0
87: PetscScalar prp,uru,vrv,wrw
88: PetscScalar p,u,v,w,ubar,c2,c
89: PetscScalar eig1,eig2,eig3,eig4,eigeps
90: PetscScalar phi1,phi2,phi3,phi4,phi5
91: PetscScalar phi6,phi7,phi8,phi9
92: PetscScalar rhs1,rhs1p,rhs1u,rhs1v,rhs1w
93: PetscScalar rhs2,rhs2p,rhs2u,rhs2v,rhs2w
94: PetscScalar rhs3,rhs3p,rhs3u,rhs3v,rhs3w
95: PetscScalar rhs4,rhs4p,rhs4u,rhs4v,rhs4w
96: PetscScalar c2inv
97: PetscScalar y11,y21,y31,y41,y12,y22,y32,y42
98: PetscScalar y13,y23,y33,y43,y14,y24,y34,y44
99: PetscScalar t11,t21,t31,t41,t12,t22,t32,t42
100: PetscScalar t13,t23,t33,t43,t14,t24,t34,t44
101: PetscScalar ti11,ti21,ti31,ti41
102: PetscScalar ti12,ti22,ti32,ti42
103: PetscScalar ti13,ti23,ti33,ti43
104: PetscScalar ti14,ti24,ti34,ti44
105: PetscScalar a11,a21,a31,a41,a12,a22,a32,a42
106: PetscScalar a13,a23,a33,a43,a14,a24,a34,a44
107: PetscScalar pb,pbp,pbu,pbv,pbw
108: PetscScalar ub,ubp,ubu,ubv,ubw
109: PetscScalar vb,vbp,vbu,vbv,vbw
110: PetscScalar wb,wbp,wbu,wbv,wbw
111: PetscScalar unormb,unormbp,unormbu
112: PetscScalar unormbv,unormbw
115: !
116: ! Loop over the edges and fill the matrix
117: ! First just zero it out
118: !
119: flops = 0.0
120: call MatZeroEntries(A,ierr)
121: ! write(6,555)res0,resc,cfl,cfl1,cfl2
122: ! 555 format(1h ,'In FILLA res0 resc cfl cfl1 cfl2 =',5(e14.7,1x))
123: !
124: #if defined(INTERLACING)
125: do i = 1,nnodes
126: temp = vol(i)/(cfl*cdt(i))
127: do j = 1,4
128: ir = j - 1 + 4*(i-1)
129: #if defined(FASTMATSET)
130: call MatSetValues4(A,1,ir,1,ir,temp)
131: #else
132: call MatSetValuesLocal(A,1,ir,1,ir,temp,ADD_VALUES,ierr)
133: #endif
134: enddo
135: enddo
136: flops = flops + 2.0*nnodes
137: #else
138: do j = 1,4
139: do i = 1,nnodes
140: temp = vol(i)/(cfl*cdt(i))
141: ir = i - 1 + nnodes*(j-1)
142: call MatSetValues(A,1,ir,1,ir,temp,ADD_VALUES,ierr)
143: enddo
144: enddo
145: flops = flops + 4.0*2*nnodes
146: #endif
148: ! call MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY,ierr)
149: ! call MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY,ierr)
150: ! print *, "Finished doing time stepping part to diagonal term"
151: !
152: ! Now fill A from interior contributions
153: ! We will fix the boundaries later
154: !
155: do 1040 n = 1, nedge
156: node1 = eptr(n,1)
157: node2 = eptr(n,2)
158: if ((node1 .le. nnodes).or.(node2 .le. nnodes)) then
159: !
160: ! Calculate unit normal to face and length of face
161: !
162: xnorm = xn(n)
163: ynorm = yn(n)
164: znorm = zn(n)
165: rlen = rl(n)
166: !
167: ! Now lets get our other 2 vectors
168: ! For first vector, use {1,0,0} and subtract off the component
169: ! in the direction of the face normal. If the inner product of
170: ! {1,0,0} is close to unity, use {0,1,0}
171: !
172: dot = xnorm
173: if (abs(dot).lt.0.95d0)then
174: X1 = 1.d0 - dot*xnorm
175: Y1 = - dot*ynorm
176: Z1 = - dot*znorm
177: else
178: dot = ynorm
179: X1 = - dot*xnorm
180: Y1 = 1.d0 - dot*ynorm
181: Z1 = - dot*znorm
182: end if
183: !
184: ! Normalize the first vector
185: !
186: size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
187: X1 = X1/size
188: Y1 = Y1/size
189: Z1 = Z1/size
190: !
191: ! Take cross-product of normal and V1 to get V2
192: !
193: X2 = ynorm*Z1 - znorm*Y1
194: Y2 = znorm*X1 - xnorm*Z1
195: Z2 = xnorm*Y1 - ynorm*X1
196: !
197: ! Variables on left
198: !
199: pL = qnode(1,node1)
200: uL = qnode(2,node1)
201: vL = qnode(3,node1)
202: wL = qnode(4,node1)
203: ubarL = xnorm*uL + ynorm*vL + znorm*wL
204: c2L = ubarL*ubarL + beta
205: cL = sqrt(c2L)
206: !
207: ! Variables on right
208: !
209: pR = qnode(1,node2)
210: uR = qnode(2,node2)
211: vR = qnode(3,node2)
212: wR = qnode(4,node2)
213: ubarR = xnorm*uR + ynorm*vR + znorm*wR
214: c2R = ubarR*ubarR + beta
215: cR = sqrt(c2R)
216: !
217: ! Regular Jacobians on left
218: !
219: dfp(1,1) = 0.0d0
220: dfp(1,2) = rlen*beta*xnorm
221: dfp(1,3) = rlen*beta*ynorm
222: dfp(1,4) = rlen*beta*znorm
223: !
224: dfp(2,1) = rlen*xnorm
225: dfp(2,2) = rlen*(ubarL + xnorm*uL)
226: dfp(2,3) = rlen*ynorm*uL
227: dfp(2,4) = rlen*znorm*uL
228: !
229: dfp(3,1) = rlen*ynorm
230: dfp(3,2) = rlen*xnorm*vL
231: dfp(3,3) = rlen*(ubarL + ynorm*vL)
232: dfp(3,4) = rlen*znorm*vL
233: !
234: dfp(4,1) = rlen*znorm
235: dfp(4,2) = rlen*xnorm*wL
236: dfp(4,3) = rlen*ynorm*wL
237: dfp(4,4) = rlen*(ubarL + znorm*wL)
238: !
239: ! Now compute eigenvalues and |A| from averaged variables
240: ! Avergage variables
241: !
242: p = .5d0*(pL + pR)
243: u = .5d0*(uL + uR)
244: v = .5d0*(vL + vR)
245: w = .5d0*(wL + wR)
246: ubar = xnorm*u + ynorm*v + znorm*w
247: c2 = ubar*ubar + beta
248: c = sqrt(c2)
249: !
250: eig1 = abs(ubar)
251: eig2 = abs(ubar)
252: eig3 = abs(ubar + c)
253: eig4 = abs(ubar - c)
254: !
255: ! Put in the eigenvalue smoothing stuff
256: !
257: eigeps = .1d0*(abs(ubar) + abs(c))
258: !
259: ! if (eig1.lt.eigeps)eig1 = .5*(eig1**2/eigeps + eigeps)
260: ! if (eig2.lt.eigeps)eig2 = .5*(eig2**2/eigeps + eigeps)
261: ! if (eig3.lt.eigeps)eig3 = .5*(eig3**2/eigeps + eigeps)
262: ! if (eig4.lt.eigeps)eig4 = .5*(eig4**2/eigeps + eigeps)
263: !
264: eig1 = rlen*eig1
265: eig2 = rlen*eig2
266: eig3 = rlen*eig3
267: eig4 = rlen*eig4
268: !
269: phi1 = xnorm*beta + u*ubar
270: phi2 = ynorm*beta + v*ubar
271: phi3 = znorm*beta + w*ubar
272: phi4 = Y2*phi3 - Z2*phi2
273: phi5 = Z2*phi1 - X2*phi3
274: phi6 = X2*phi2 - Y2*phi1
275: phi7 = Z1*phi2 - Y1*phi3
276: phi8 = X1*phi3 - Z1*phi1
277: phi9 = Y1*phi1 - X1*phi2
278: !
279: ! Components of T(inverse) (call this y)
280: !
281: c2inv = 1.d0/c2
282: y11 = -c2inv*(u*phi4 + v*phi5 + w*phi6)/beta
283: y21 = -c2inv*(u*phi7 + v*phi8 + w*phi9)/beta
284: y31 = .5d0*c2inv*(c-ubar)/beta
285: y41 = -.5d0*c2inv*(c+ubar)/beta
287: y12 = c2inv*phi4
288: y22 = c2inv*phi7
289: y32 = c2inv*.5d0*xnorm
290: y42 = c2inv*.5d0*xnorm
292: y13 = c2inv*phi5
293: y23 = c2inv*phi8
294: y33 = c2inv*.5d0*ynorm
295: y43 = c2inv*.5d0*ynorm
297: y14 = c2inv*phi6
298: y24 = c2inv*phi9
299: y34 = c2inv*.5d0*znorm
300: y44 = c2inv*.5d0*znorm
301: !
302: ! Now get elements of T
303: !
304: t11 = 0.0d0
305: t21 = X1
306: t31 = Y1
307: t41 = Z1
309: t12 = 0.0d0
310: t22 = X2
311: t32 = Y2
312: t42 = Z2
314: t13 = c*beta
315: t23 = xnorm*beta + u*(ubar + c)
316: t33 = ynorm*beta + v*(ubar + c)
317: t43 = znorm*beta + w*(ubar + c)
319: t14 = -c*beta
320: t24 = xnorm*beta + u*(ubar - c)
321: t34 = ynorm*beta + v*(ubar - c)
322: t44 = znorm*beta + w*(ubar - c)
323: !
324: ! Compute T*|lambda|*T(inv)
325: !
326: a11 = eig1*t11*y11 + eig2*t12*y21 + eig3*t13*y31 + eig4*t14*y41
327: a12 = eig1*t11*y12 + eig2*t12*y22 + eig3*t13*y32 + eig4*t14*y42
328: a13 = eig1*t11*y13 + eig2*t12*y23 + eig3*t13*y33 + eig4*t14*y43
329: a14 = eig1*t11*y14 + eig2*t12*y24 + eig3*t13*y34 + eig4*t14*y44
331: a21 = eig1*t21*y11 + eig2*t22*y21 + eig3*t23*y31 + eig4*t24*y41
332: a22 = eig1*t21*y12 + eig2*t22*y22 + eig3*t23*y32 + eig4*t24*y42
333: a23 = eig1*t21*y13 + eig2*t22*y23 + eig3*t23*y33 + eig4*t24*y43
334: a24 = eig1*t21*y14 + eig2*t22*y24 + eig3*t23*y34 + eig4*t24*y44
336: a31 = eig1*t31*y11 + eig2*t32*y21 + eig3*t33*y31 + eig4*t34*y41
337: a32 = eig1*t31*y12 + eig2*t32*y22 + eig3*t33*y32 + eig4*t34*y42
338: a33 = eig1*t31*y13 + eig2*t32*y23 + eig3*t33*y33 + eig4*t34*y43
339: a34 = eig1*t31*y14 + eig2*t32*y24 + eig3*t33*y34 + eig4*t34*y44
341: a41 = eig1*t41*y11 + eig2*t42*y21 + eig3*t43*y31 + eig4*t44*y41
342: a42 = eig1*t41*y12 + eig2*t42*y22 + eig3*t43*y32 + eig4*t44*y42
343: a43 = eig1*t41*y13 + eig2*t42*y23 + eig3*t43*y33 + eig4*t44*y43
344: a44 = eig1*t41*y14 + eig2*t42*y24 + eig3*t43*y34 + eig4*t44*y44
345: !
346: ! Form .5*(A + |A|)
347: !
348: dfp(1,1) = .5d0*(dfp(1,1) + a11)
349: dfp(1,2) = .5d0*(dfp(1,2) + a12)
350: dfp(1,3) = .5d0*(dfp(1,3) + a13)
351: dfp(1,4) = .5d0*(dfp(1,4) + a14)
352: !
353: dfp(2,1) = .5d0*(dfp(2,1) + a21)
354: dfp(2,2) = .5d0*(dfp(2,2) + a22)
355: dfp(2,3) = .5d0*(dfp(2,3) + a23)
356: dfp(2,4) = .5d0*(dfp(2,4) + a24)
357: !
358: dfp(3,1) = .5d0*(dfp(3,1) + a31)
359: dfp(3,2) = .5d0*(dfp(3,2) + a32)
360: dfp(3,3) = .5d0*(dfp(3,3) + a33)
361: dfp(3,4) = .5d0*(dfp(3,4) + a34)
362: !
363: dfp(4,1) = .5d0*(dfp(4,1) + a41)
364: dfp(4,2) = .5d0*(dfp(4,2) + a42)
365: dfp(4,3) = .5d0*(dfp(4,3) + a43)
366: dfp(4,4) = .5d0*(dfp(4,4) + a44)
367: !
368: ! Now get regular Jacobians on right
369: !
370: dfm(1,1) = 0.0d0
371: dfm(1,2) = rlen*beta*xnorm
372: dfm(1,3) = rlen*beta*ynorm
373: dfm(1,4) = rlen*beta*znorm
374: !
375: dfm(2,1) = rlen*xnorm
376: dfm(2,2) = rlen*(ubarR + xnorm*uR)
377: dfm(2,3) = rlen*ynorm*uR
378: dfm(2,4) = rlen*znorm*uR
379: !
380: dfm(3,1) = rlen*ynorm
381: dfm(3,2) = rlen*xnorm*vR
382: dfm(3,3) = rlen*(ubarR + ynorm*vR)
383: dfm(3,4) = rlen*znorm*vR
384: !
385: dfm(4,1) = rlen*znorm
386: dfm(4,2) = rlen*xnorm*wR
387: dfm(4,3) = rlen*ynorm*wR
388: dfm(4,4) = rlen*(ubarR + znorm*wR)
389: !
390: ! Form .5*(A - |A|)
391: !
392: dfm(1,1) = .5d0*(dfm(1,1) - a11)
393: dfm(1,2) = .5d0*(dfm(1,2) - a12)
394: dfm(1,3) = .5d0*(dfm(1,3) - a13)
395: dfm(1,4) = .5d0*(dfm(1,4) - a14)
396: !
397: dfm(2,1) = .5d0*(dfm(2,1) - a21)
398: dfm(2,2) = .5d0*(dfm(2,2) - a22)
399: dfm(2,3) = .5d0*(dfm(2,3) - a23)
400: dfm(2,4) = .5d0*(dfm(2,4) - a24)
401: !
402: dfm(3,1) = .5d0*(dfm(3,1) - a31)
403: dfm(3,2) = .5d0*(dfm(3,2) - a32)
404: dfm(3,3) = .5d0*(dfm(3,3) - a33)
405: dfm(3,4) = .5d0*(dfm(3,4) - a34)
406: !
407: dfm(4,1) = .5d0*(dfm(4,1) - a41)
408: dfm(4,2) = .5d0*(dfm(4,2) - a42)
409: dfm(4,3) = .5d0*(dfm(4,3) - a43)
410: dfm(4,4) = .5d0*(dfm(4,4) - a44)
412: flops = flops + 465.0
413: !
414: ! Now take care of contribution to node 1
415: !
416: ! idiag = iau(node1)
417: !
418: ! Diagonal piece
419: !
420: if (node1 .le. nnodes) then
421: #if defined(INTERLACING)
422: #if defined(BLOCKING)
423: irow(1) = node1 - 1
424: icol(1) = node1 - 1
425: icol(2) = node2 - 1
426: #if defined(FASTMATSET)
427: call MatSetValuesBlocked4(A,1,irow,2,icol,val1)
428: #else
429: call MatSetValuesBlockedLocal(A,1,irow,2,icol, &
430: & val1,ADD_VALUES,ierr)
431: #endif
432: #else
433: do j = 1,4
434: irow(j) = 4*(node1-1)+j-1
435: icol(j) = irow(j)
436: icol(4+j) = 4*(node2-1)+j-1
437: enddo
438: call MatSetValuesLocal(A,4,irow,8,icol,val1,ADD_VALUES,ierr)
439: #endif
440: #else
441: do j = 1,4
442: irow(j) = (node1-1)+(j-1)*nnodes
443: icol(j) = irow(j)
444: icol(4+j) = (node2-1)+(j-1)*nnodes
445: enddo
446: call MatSetValues(A,4,irow,8,icol,val1,ADD_VALUES,ierr)
447: #endif
448: endif
450: !
451: ! Now do the second node
452: !
453: if (node2 .le. nnodes) then
454: !
455: ! Exchange elements in place
456: do j = 1,4
457: do k = 1,8
458: ! temp = -val1(k,j)
459: ! val1(k,j) = -val1(k+4,j)
460: ! val1(k+4,j) = temp
461: val1(k,j) = -val1(k,j)
462: enddo
463: enddo
464: !
465: ! call CHK_ERR(irank,ierr,irow(1),icol(1))
466: #if defined(INTERLACING)
467: #if defined(BLOCKING)
468: irow(1) = node2 - 1
469: icol(1) = node1 - 1
470: icol(2) = node2 - 1
471: #if defined(FASTMATSET)
472: call MatSetValuesBlocked4(A,1,irow,2,icol,val1)
473: #else
474: call MatSetValuesBlockedLocal(A,1,irow,2,icol, &
475: & val1,ADD_VALUES,ierr)
476: #endif
477: #else
478: do j = 1,4
479: irow(j) = 4*(node2-1)+j-1
480: icol(j) = 4*(node1-1)+j-1
481: icol(4+j) = irow(j)
482: enddo
483: call MatSetValuesLocal(A,4,irow,8,icol,val1,ADD_VALUES,ierr)
484: #endif
485: #else
486: do j = 1,4
487: irow(j) = (node2-1)+(j-1)*nnodes
488: icol(j) = (node1-1)+(j-1)*nnodes
489: icol(4+j) = irow(j)
490: enddo
491: call MatSetValues(A,4,irow,8,icol,val1,ADD_VALUES,ierr)
492: #endif
493: endif
495: endif
496: 1040 continue
497: !
498: ! Now loop over the boundaries
499: !
500: ! For inviscid surface add contribution from pressure
501: !
502: do 1070 i = 1,nsnode
503: inode = isnode(i)
504: if (inode .le. nnodes) then
505: xnorm = sxn(i)
506: ynorm = syn(i)
507: znorm = szn(i)
508: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
509: xnorm = xnorm/area
510: ynorm = ynorm/area
511: znorm = znorm/area
512: !
513: val(1) = area*xnorm
514: val(2) = area*ynorm
515: val(3) = area*znorm
516: #if defined(INTERLACING)
517: irow(1) = 4*(inode-1) + 1
518: irow(2) = 4*(inode-1) + 2
519: irow(3) = 4*(inode-1) + 3
520: ic = 4*(inode - 1)
521: #if defined(FASTMATSET)
522: call MatSetValues4(A,3,irow,1,ic,val)
523: #else
524: call MatSetValuesLocal(A,3,irow,1,ic,val,ADD_VALUES,ierr)
525: #endif
526: #else
527: irow(1) = inode - 1 + nnodes
528: irow(2) = inode - 1 + nnodes*2
529: irow(3) = inode - 1 + nnodes*3
530: ic = inode - 1
531: call MatSetValues(A,3,irow,1,ic,val,ADD_VALUES,ierr)
532: #endif
533: flops = flops + 12.0
534: endif
536: !
537: ! idiag = iau(inode)
538: ! A(idiag,2,1) = A(idiag,2,1) + area*xnorm
539: ! A(idiag,3,1) = A(idiag,3,1) + area*ynorm
540: ! A(idiag,4,1) = A(idiag,4,1) + area*znorm
541: 1070 continue
542: ! print *, "Finished doing inviscid nodes"
543: !
544: ! Now do viscous faces
545: !
546: ! do 1080 i = 1,nvnode
547: ! inode = ivnode(i)
548: ! idiag = iau(inode)
549: !
550: ! First zero out all the ones on the row and then
551: ! refill them so that the velocity is just zero on body
552: !
553: ! jstart = ia(inode)
554: ! jend = ia(inode+1) - 1
555: !
556: ! do 1060 j=jstart,jend
557: !
558: ! If this is not a diagonal zero it out
559: ! (This way we dont disturb the row for the coninuity equation
560: !
561: ! if (j.ne.idiag)then
562: ! A(j,1,1) = 0.0
563: ! A(j,1,2) = 0.0
564: ! A(j,1,3) = 0.0
565: ! A(j,1,4) = 0.0
566: !
567: ! A(j,2,1) = 0.0
568: ! A(j,2,2) = 0.0
569: ! A(j,2,3) = 0.0
570: ! A(j,2,4) = 0.0
571: !
572: ! A(j,3,1) = 0.0
573: ! A(j,3,2) = 0.0
574: ! A(j,3,3) = 0.0
575: ! A(j,3,4) = 0.0
576: !
577: ! A(j,4,1) = 0.0
578: ! A(j,4,2) = 0.0
579: ! A(j,4,3) = 0.0
580: ! A(j,4,4) = 0.0
581: !
582: ! end if
583: !1060 continue
584: !
585: ! Now set the diagonal for the momentum equations
586: !
587: ! A(idiag,2,1) = 0.0
588: ! A(idiag,2,2) = 1.0
589: ! A(idiag,2,3) = 0.0
590: ! A(idiag,2,4) = 0.0
591: !
592: ! A(idiag,3,1) = 0.0
593: ! A(idiag,3,2) = 0.0
594: ! A(idiag,3,3) = 1.0
595: ! A(idiag,3,4) = 0.0
596: !
597: ! A(idiag,4,1) = 0.0
598: ! A(idiag,4,2) = 0.0
599: ! A(idiag,4,3) = 0.0
600: ! A(idiag,4,4) = 1.0
601: !
602: !1080 continue
603: !
604: ! Far-field boundary points
605: !
606: do 1090 i = 1,nfnode
607: inode = ifnode(i)
608: if (inode .le. nnodes) then
609: xnorm = fxn(i)
610: ynorm = fyn(i)
611: znorm = fzn(i)
612: area = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
613: xnorm = xnorm/area
614: ynorm = ynorm/area
615: znorm = znorm/area
616: !
617: ! Now lets get our other 2 vectors
618: ! For first vector, use {1,0,0} and subtract off the component
619: ! in the direction of the face normal. If the inner product of
620: ! {1,0,0} is close to unity, use {0,1,0}
621: !
622: dot = xnorm
623: if (abs(dot).lt.0.95d0)then
624: X1 = 1.d0 - dot*xnorm
625: Y1 = - dot*ynorm
626: Z1 = - dot*znorm
627: else
628: dot = ynorm
629: X1 = - dot*xnorm
630: Y1 = 1.d0 - dot*ynorm
631: Z1 = - dot*znorm
632: end if
633: !
634: ! Normalize the first vector (V1)
635: !
636: size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
637: X1 = X1/size
638: Y1 = Y1/size
639: Z1 = Z1/size
640: !
641: ! Take cross-product of normal with V1 to get V2
642: !
643: X2 = ynorm*Z1 - znorm*Y1
644: Y2 = znorm*X1 - xnorm*Z1
645: Z2 = xnorm*Y1 - ynorm*X1
646: !
647: ! Calculate elements of T and T(inverse) evaluated at freestream
648: !
649: ubar0 = xnorm*u0 + ynorm*v0 + znorm*w0
650: c20 = ubar0*ubar0 + beta
651: c0 = sqrt(c20)
652: phi1 = xnorm*beta + u0*ubar0
653: phi2 = ynorm*beta + v0*ubar0
654: phi3 = znorm*beta + w0*ubar0
655: phi4 = Y2*phi3 - Z2*phi2
656: phi5 = Z2*phi1 - X2*phi3
657: phi6 = X2*phi2 - Y2*phi1
658: phi7 = Z1*phi2 - Y1*phi3
659: phi8 = X1*phi3 - Z1*phi1
660: phi9 = Y1*phi1 - X1*phi2
662: t11 = 0.0d0
663: t21 = X1
664: t31 = Y1
665: t41 = Z1
667: t12 = 0.0d0
668: t22 = X2
669: t32 = Y2
670: t42 = Z2
672: t13 = c0*beta
673: t23 = xnorm*beta + u0*(ubar0 + c0)
674: t33 = ynorm*beta + v0*(ubar0 + c0)
675: t43 = znorm*beta + w0*(ubar0 + c0)
677: t14 = -c0*beta
678: t24 = xnorm*beta + u0*(ubar0 - c0)
679: t34 = ynorm*beta + v0*(ubar0 - c0)
680: t44 = znorm*beta + w0*(ubar0 - c0)
682: ti11 = -(u0*phi4 + v0*phi5 + w0*phi6)/beta/c20
683: ti21 = -(u0*phi7 + v0*phi8 + w0*phi9)/beta/c20
684: ti31 = (c0 - ubar0)/(2.d0*beta*c20)
685: ti41 = -(c0 + ubar0)/(2.d0*beta*c20)
687: ti12 = phi4/c20
688: ti22 = phi7/c20
689: ti32 = .5d0*xnorm/c20
690: ti42 = .5d0*xnorm/c20
692: ti13 = phi5/c20
693: ti23 = phi8/c20
694: ti33 = .5d0*ynorm/c20
695: ti43 = .5d0*ynorm/c20
697: ti14 = phi6/c20
698: ti24 = phi9/c20
699: ti34 = .5d0*znorm/c20
700: ti44 = .5d0*znorm/c20
701: !
702: ! Now, get the variables on the "inside"
703: !
704: pi = qnode(1,inode)
705: ui = qnode(2,inode)
706: vi = qnode(3,inode)
707: wi = qnode(4,inode)
708: unorm = xnorm*ui + ynorm*vi + znorm*wi
709: !
710: ! If ubar is negative, take the reference condition from outside
711: !
712: !
713: if (unorm.gt.0.0d0)then
714: pr = pi
715: prp = 1.0d0
716: ur = ui
717: uru = 1.0d0
718: vr = vi
719: vrv = 1.0d0
720: wr = wi
721: wrw = 1.0d0
722: else
723: pr = p0
724: prp = 0.0d0
725: ur = u0
726: uru = 0.0d0
727: vr = v0
728: vrv = 0.0d0
729: wr = w0
730: wrw = 0.0d0
731: end if
732: !
733: ! Set rhs
734: !
735: rhs1 = ti11*pr + ti12*ur + ti13*vr + ti14*wr
736: rhs1p = ti11*prp
737: rhs1u = ti12*uru
738: rhs1v = ti13*vrv
739: rhs1w = ti14*wrw
740: rhs2 = ti21*pr + ti22*ur + ti23*vr + ti24*wr
741: rhs2p = ti21*prp
742: rhs2u = ti22*uru
743: rhs2v = ti23*vrv
744: rhs2w = ti24*wrw
745: rhs3 = ti31*pi + ti32*ui + ti33*vi + ti34*wi
746: rhs3p = ti31
747: rhs3u = ti32
748: rhs3v = ti33
749: rhs3w = ti34
750: rhs4 = ti41*p0 + ti42*u0 + ti43*v0 + ti44*w0
751: rhs4p = 0.0d0
752: rhs4u = 0.0d0
753: rhs4v = 0.0d0
754: rhs4w = 0.0d0
755: !
756: ! Now do matrix multiplication to get values on boundary
757: !
758: pb = t13*rhs3 + t14*rhs4
759: pbp = t13*rhs3p !+ t14*rhs4p
760: pbu = t13*rhs3u !+ t14*rhs4u
761: pbv = t13*rhs3v !+ t14*rhs4v
762: pbw = t13*rhs3w !+ t14*rhs4w
763: ub = t21*rhs1 + t22*rhs2 + t23*rhs3 + t24*rhs4
764: ubp = t21*rhs1p + t22*rhs2p + t23*rhs3p !+ t24*rhs4p
765: ubu = t21*rhs1u + t22*rhs2u + t23*rhs3u !+ t24*rhs4u
766: ubv = t21*rhs1v + t22*rhs2v + t23*rhs3v !+ t24*rhs4v
767: ubw = t21*rhs1w + t22*rhs2w + t23*rhs3w !+ t24*rhs4w
768: vb = t31*rhs1 + t32*rhs2 + t33*rhs3 + t34*rhs4
769: vbp = t31*rhs1p + t32*rhs2p + t33*rhs3p !+ t34*rhs4p
770: vbu = t31*rhs1u + t32*rhs2u + t33*rhs3u !+ t34*rhs4u
771: vbv = t31*rhs1v + t32*rhs2v + t33*rhs3v !+ t34*rhs4v
772: vbw = t31*rhs1w + t32*rhs2w + t33*rhs3w !+ t34*rhs4w
773: wb = t41*rhs1 + t42*rhs2 + t43*rhs3 + t44*rhs4
774: wbp = t41*rhs1p + t42*rhs2p + t43*rhs3p !+ t44*rhs4p
775: wbu = t41*rhs1u + t42*rhs2u + t43*rhs3u !+ t44*rhs4u
776: wbv = t41*rhs1v + t42*rhs2v + t43*rhs3v !+ t44*rhs4v
777: wbw = t41*rhs1w + t42*rhs2w + t43*rhs3w !+ t44*rhs4w
779: unormb = xnorm*ub + ynorm*vb + znorm*wb
780: unormbp = xnorm*ubp + ynorm*vbp + znorm*wbp
781: unormbu = xnorm*ubu + ynorm*vbu + znorm*wbu
782: unormbv = xnorm*ubv + ynorm*vbv + znorm*wbv
783: unormbw = xnorm*ubw + ynorm*vbw + znorm*wbw
785: !
786: ! Now add contribution to lhs
787: !
789: val(1) = area*beta*unormbp
790: val(2) = area*beta*unormbu
791: val(3) = area*beta*unormbv
792: val(4) = area*beta*unormbw
793: !
794: val(5) = area*(ub*unormbp + unormb*ubp + xnorm*pbp)
795: val(6) = area*(ub*unormbu + unormb*ubu + xnorm*pbu)
796: val(7) = area*(ub*unormbv + unormb*ubv + xnorm*pbv)
797: val(8) = area*(ub*unormbw + unormb*ubw + xnorm*pbw)
798: !
799: val(9) = area*(vb*unormbp + unormb*vbp + ynorm*pbp)
800: val(10) = area*(vb*unormbu + unormb*vbu + ynorm*pbu)
801: val(11) = area*(vb*unormbv + unormb*vbv + ynorm*pbv)
802: val(12) = area*(vb*unormbw + unormb*vbw + ynorm*pbw)
803: !
804: val(13) = area*(wb*unormbp + unormb*wbp + znorm*pbp)
805: val(14) = area*(wb*unormbu + unormb*wbu + znorm*pbu)
806: val(15) = area*(wb*unormbv + unormb*wbv + znorm*pbv)
807: val(16) = area*(wb*unormbw + unormb*wbw + znorm*pbw)
808: !
809: #if defined(INTERLACING)
810: #if defined(BLOCKING)
811: irow(1) = inode - 1
812: #if defined(FASTMATSET)
813: call MatSetValuesBlocked4(A,1,irow,1,irow,val)
814: #else
815: call MatSetValuesBlockedLocal(A,1,irow,1,irow,val, &
816: & ADD_VALUES,ierr)
817: #endif
818: #else
819: do k = 1,4
820: irow(k) = 4*(inode-1)+k-1
821: enddo
822: call MatSetValuesLocal(A,4,irow,4,irow,val,ADD_VALUES,ierr)
823: #endif
824: #else
825: do k = 1,4
826: irow(k) = inode - 1 + nnodes*(k-1)
827: enddo
828: call MatSetValues(A,4,irow,4,irow,val,ADD_VALUES,ierr)
829: #endif
831: flops = flops + 337.0
832: endif
833: !
834: 1090 continue
836: ! print *, "Finished doing far field nodes"
838: ! Assemble matrix
839: call PetscLogFlops(flops,ierr)
841: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
842: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
843: ! call MatView(A, PETSC_VIEWER_STDOUT_SELF,ierr)
844: flag = SAME_NONZERO_PATTERN
845: !
846: ! End of subroutine FILLA
847: !
848: return
849: end
850: !
851: !
853: subroutine CHK_ERR(irank, ierr, irow, icol)
854: implicit none
855: #include <petsc/finclude/petscsys.h>
856: integer irank,ierr,irow,icol
857: if (ierr .gt. 0) then
858: write(*,*) 'On processor ',irank, ': Non-zero entry in row ', &
859: & irow, ' and column ',icol,' is beyond the pre-allocated memory'
860: endif
861: return
862: end