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