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