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