Actual source code: user.F

      1: !---------------------------------------------------------------
      2: ! The following subroutines are from node2t.f in the original
      3: ! code - D. K. Kaushik (1/17/97)
      4: !---------------------------------------------------------------
      5: !
      6: !
      7: ! 2-D Navier Stokes on Unstructured Grids

      9: !============================== FORLINK ==============================72
     10: !
     11: !  FORLINK establishes links between FORTRAN common blocks and C
     12: !
     13: !=====================================================================72
     14:       subroutine FORLINK()

     16:       implicit none
     17: #include <petsc/finclude/petscsys.h>
     18:       PetscScalar  title(20),beta,alpha
     19:       PetscScalar  Re,dt,tot,res0,resc
     20:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
     21:       PetscScalar  cfl1,cfl2
     22:       integer nsmoth,iflim,itran,nbtran,jupdate,                        &
     23:      &        nstage,ncyct,iramp,nitfo
     24:       PetscScalar  gtol
     25:       integer icycle,nsrch,ilu0,ifcn
     26:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                 &
     27:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
     28:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,         &
     29:      &             nstage,ncyct,iramp,nitfo
     30:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn

     32:       call CLINK(title,cfl1,gtol)
     33: !
     34: ! End of subroutine FORLINK
     35: !
     36:       return
     37:       end

     39: !============================== Block_Initialization  ======== =======72
     40: !
     41: ! Initializes the common blocks members for turbulence model
     42: !
     43: !=====================================================================72
     44:       block data Block_Initialization
     45:       implicit none
     46: #include <petsc/finclude/petscsys.h>
     47:       PetscScalar vkar,cmu,ce1,ce2
     48:       PetscScalar aplus1,aplus2,turbinf
     49:       PetscScalar cb1,sig,cb2,cw1,cw2
     50:       PetscScalar cw3,cv1,ct1,ct2,ct3,ct4
     51:       common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
     52:       common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
     53:       data vkar,cmu,ce1,ce2/0.41,0.09,1.2,2.0 /
     54:       data aplus1,aplus2,turbinf/26.0,10.0,0.1 /
     55:       data cb1,sig,cb2,cw2,cw3/0.1355,0.66667,0.622,0.3,2.0/
     56: !
     57: ! Comment out old coefficients
     58: !      data cv1,ct1,ct2,ct3,ct4/7.1,1.0,2.0,1.1,2.0/
     59: !
     60:        data cv1,ct1,ct2,ct3,ct4/7.1d0,1.0d0,2.0d0,1.2d0,0.5d0/

     62:       end

     64: !================================== INIT =============================72
     65: !
     66: ! Initializes the flow field
     67: !
     68: !=====================================================================72
     69:       subroutine INIT(nnodes, qvec, turbre, amut,nvnode, ivnode,irank)
     70:       implicit none
     71: #include <petsc/finclude/petscsys.h>
     72:       PetscScalar  turbre(1),amut(1)
     73:       integer ivnode(1),nnodes,nvnode,irank
     74:       PetscScalar  title(20),beta,alpha
     75:       PetscScalar  Re,dt,tot,res0,resc
     76:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
     77:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
     78:       PetscScalar  p0,rho0,c0,u0,v0,w0,et0,h0,pt0
     79:       PetscScalar  vkar,cmu,ce1,ce2
     80:       PetscScalar  aplus1,aplus2,turbinf
     81:       PetscScalar  cb1,sig,cb2,cw1,cw2,cw3
     82:       PetscScalar  cv1,ct1,ct2,ct3,ct4
     83:       PetscScalar  cfl1,cfl2
     84:       integer nsmoth,iflim,itran,nbtran,jupdate,                          &
     85:      &        nstage,ncyct,iramp,nitfo
     86:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                   &
     87:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
     88:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
     89:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
     90:       common/turb/vkar,cmu,ce1,ce2,aplus1,aplus2,turbinf
     91:       common/spalrt/cb1,sig,cb2,cw1,cw2,cw3,cv1,ct1,ct2,ct3,ct4
     92:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,           &
     93:      &             nstage,ncyct,iramp,nitfo
     94:       PetscScalar  pi,conv,rmu,chi,fv1
     95:       Integer i,k,n
     96: #if defined(INTERLACING)
     97:        PetscScalar qvec(4,nnodes)
     98: #define qnod(i,j) qvec(i,j)
     99: #else
    100:        PetscScalar qvec(nnodes,4)
    101: #define qnod(i,j) qvec(j,i)
    102: #endif
    103: !
    104:       cw1 = cb1/vkar/vkar + (1.0d0 + cb2)/sig
    105:       if (ivisc.gt.4d0)turbinf = 1.341946d0
    106: !      if (ivisc.gt.4)turbinf = 0.5
    107:       if (ivisc.gt.4.and.itran.eq.1)turbinf = 0.1d0
    108: !
    109: ! Note that for Spalarts model, I use turbinf as the freestream value of
    110: ! the dependent variable just as in the Baldwin-Barth model.
    111: ! The constant is set so that in the freestream, nu_t=0.009 (1.341946)
    112: !
    113: !     print *, "I am in INIT"
    114:       res0 = 1.0d0
    115:       resc = 1.0d0

    117:       gamma = 1.0d0
    118:       gm1   = 1.0d0
    119:       gp1   = 1.0d0
    120:       gm1g  = 1.0d0
    121:       gp1g  = 1.0d0
    122:       ggm1  = 1.0d0

    124:       pi = 4.0d0*datan(1.0d0)
    125:       conv = 180.0d0/pi
    126:       p0    = 1.0d0
    127: #if defined(CFL3D_AXIS)
    128:       u0    = cos(alpha/conv)
    129:       v0    = 0.0d0
    130:       w0    = sin(alpha/conv)
    131: #else
    132:       u0    = cos(alpha/conv)
    133:       v0    = sin(alpha/conv)
    134:       w0    = 0.0d0
    135: #endif
    136:       do n = 1,nnodes
    137:          qnod(1,n) = p0
    138:          qnod(2,n) = u0
    139:          qnod(3,n) = v0
    140:          qnod(4,n) = w0
    141:       enddo
    142:       if (ivisc.eq.3) then
    143:        do n = 1,nnodes
    144:          turbre(n)  = 0.0d0
    145:          amut(n)    = 0.0d0
    146:        enddo
    147:       endif
    148: !
    149: ! If viscous, zero out the velocity on the surface
    150: !
    151: !     print *, "Just Before Viscous"

    153:       do 9010 i = 1,nvnode
    154: !
    155: ! Compute the velocity normal to the surface
    156: !
    157:         k       = ivnode(i)
    158:         qnod(2,k) = 0.0d0
    159:         qnod(3,k) = 0.0d0
    160:         qnod(4,k) = 0.0d0
    161:  9010 continue
    162: !
    163: ! If turbulent, initialize turbre
    164: !
    165:       if (ivisc.eq.3.or.ivisc.eq.4) then
    166:          if (irank .eq. 0) write(10,110)
    167:          if (irank .eq. 0) write(10,120)vkar,cmu,ce1,ce2
    168:          if (irank .eq. 0) write(10,130)aplus1,aplus2,turbinf
    169:          do 1010 n = 1,nnodes
    170:             turbre(n) = turbinf
    171:             amut(n)   = cmu*turbinf
    172:  1010    continue
    173:       end if

    175:       if (ivisc.eq.5.or.ivisc.eq.6) then
    176:          if (irank .eq. 0) write(10,110)
    177:          if (irank .eq. 0) write(10,140)vkar,cb1,sig,cb2
    178:          if (irank .eq. 0) write(10,150)cw1,cw2,cw3,cv1
    179:          if (irank .eq. 0) write(10,160)ct1,ct2,ct3,ct4
    180:          do 1020 n = 1,nnodes
    181:             turbre(n) = turbinf
    182:             rmu = 1.0d0
    183:             chi = turbre(n)/rmu
    184:             fv1 = chi**3/(chi**3 + cv1**3)
    185:             amut(n)   = fv1*turbinf
    186:  1020    continue
    187:       end if
    188: !     print *, "I am out of INIT"

    190:       return
    191:   110 format(1h ,'Parameters for turbulence model')
    192:   120 format(1h ,'k=',f10.5,' cmu=',f10.5,' ce1=',f10.5,'ce2=',f10.5)
    193:   130 format(1h ,'aplus1',f10.5,' aplus2=',f10.5,' turbinf=',f10.5)
    194:   140 format(1h ,'k=',f10.5,' cb1=',f10.5,' sig=',f10.5,'cb2=',f10.5)
    195:   150 format(1h ,'cw1=',f10.5,' cw2=',f10.5,' cw3=',f10.5,' cv1=',f10.5)
    196:   160 format(1h ,'ct1=',f10.5,' ct2=',f10.5,' ct3=',f10.5,' ct4=',f10.5)
    197: !
    198: ! End of subroutine INIT
    199: !

    201:       end

    204: !================================ READR1 ====================================
    205: !
    206: !  Reads input parameters
    207: !
    208: !============================================================================
    209:       subroutine READR1(ileast, irank)

    211:       implicit none
    212: #include <petsc/finclude/petscsys.h>
    213:       PetscScalar  title(20),beta,alpha
    214:       PetscScalar  Re,dt,tot,res0,resc
    215:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    216:       PetscScalar  cfl1,cfl2
    217:       integer nsmoth,iflim,itran,nbtran,jupdate,                                &
    218:      &        nstage,ncyct,iramp,nitfo
    219:       PetscScalar  gtol
    220:       integer icycle,nsrch,ilu0,ifcn
    221:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                         &
    222:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    223:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,                 &
    224:      &             nstage,ncyct,iramp,nitfo
    225:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
    226:       integer ileast, irank
    227:       integer i

    229:       read(7,10)(title(i),i=1,20)
    230:       if (irank .eq. 0) write(10,11)(title(i),i=1,20)

    232:       read(7,10)

    234:       read(7,24)mseq,ihane,ivisc,ileast,iflim,jupdate,ifcn
    235:       if (irank .eq. 0) write(10,25)mseq,ihane,ivisc
    236:       if (irank .eq. 0) write(10,28)ileast,iflim,jupdate

    238:       read(7,10)

    240:       read(7,12)beta,alpha,Re
    241:       if (irank .eq. 0) write(10,13)beta,alpha,Re

    243:       read(7,10)

    245:       read(7,26)cfl2,dt,irest,itran,nbtran
    246:       if (irank .eq. 0) write(10,27)cfl2,dt,irest

    248:       if (irank .eq. 0) then
    249:        if (ivisc.eq.5.or.ivisc.eq.6) write(10,123)itran,nbtran
    250:       endif

    252:    10 format(20a4)
    253:    11 format(1h ,20a4)
    254:    12 format(2f10.5,e14.7)
    255:    13 format(1h ,'beta = ',f10.5,' Alpha = ',f10.5,' Re = ',e14.7)
    256:    24 format(i10,i10,i10,i10,i10,i10,i10)
    257:    25 format(1h ,'mseq = ',i3,' ihane = ',i3,' ivisc=',i3)
    258:    26 format(2f10.5,3i10)
    259:    27 format(1h ,' cfl2= ',e14.7,' dt= ',f10.5,'irest= ',i5)
    260:    28 format(1h ,'ileast= ',i5,' iflim= ',i5,' jupdate= ',i5)
    261:   123 format(1h ,'itran = ',i5,' nbtran= ',i5)

    263:       return
    264:       end

    268: !================================ TECFLO =============================72
    269: !
    270: !  Writes a formatted file that contains an input file for TECPLOT
    271: !
    272: !=====================================================================72
    273:       subroutine TECFLO(nnodes,                                            &
    274:      &                  nnbound,nvbound,nfbound,                           &
    275:      &                  nnfacet,nvfacet,nffacet,                           &
    276:      &                  nsnode,nvnode,nfnode,                              &
    277:      &                  title,x,y,z,q,                                     &
    278:      &                  nnpts,nntet,nvpts,nvtet,nfpts,nftet,               &
    279:      &                  f2ntn,f2ntv,f2ntf,isnode,ivnode,ifnode,            &
    280:      &                  irank)
    281: !
    282:       implicit none
    283: #include <petsc/finclude/petscsys.h>
    284:       PetscScalar pinf
    285:       PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
    286:       PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    287:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
    288:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    289:       integer nnodes,nnbound,nvbound,nfbound,                              &
    290:      &        nnfacet,nvfacet,nffacet,                                     &
    291:      &        nsnode,nvnode,nfnode,irank
    292:       integer nntet(1),nnpts(1)
    293:       integer nvtet(1),nvpts(1)
    294:       integer nftet(1),nfpts(1)
    295:       integer f2ntn(nnfacet,4),f2ntv(nvfacet,4),f2ntf(nffacet,4)
    296:       integer isnode(1),ivnode(1),ifnode(1)
    297:       PetscScalar  x(1),y(1),z(1),q(4,nnodes)
    298:       integer i,j,i1,i2,i3,ic,isp,ist

    301:       character c4*4,title*20

    303:       pinf = p0
    304: !
    305: !   + write tecplot header
    306: !
    307:       write(13,'(a)') 'TITLE="' // title // '"'
    308:       write(13,'(a)') 'VARIABLES="X     ","Y     ","Z     ","RHO   ",',     &
    309:      &     '"U     ","V     ","W     ","P/Pinf"'
    310: !
    311: !   + do a zone-title, so we can keep track of things
    312: !   + start with solid-wall boundary surfaces
    313: !
    314:       isp   = 1
    315:       ist   = 1
    316:       do 10 i=1,nnbound
    317:         write(c4,"(i4)") i
    318:         if (i .ge.    0 .and. i .le.    9) ic = 4
    319:         if (i .ge.   10 .and. i .le.   99) ic = 3
    320:         if (i .ge.  100 .and. i .le.  999) ic = 2
    321:         if (i .ge. 1000 .and. i .le. 9999) ic = 1

    323:         write(13,1000) "nn." // c4(ic:4),nnpts(i),nntet(i)

    325:         write(13,1010) (x(isnode(j))      ,j=isp,isp+nnpts(i)-1)
    326:         write(13,1010) (y(isnode(j))      ,j=isp,isp+nnpts(i)-1)
    327:         write(13,1010) (z(isnode(j))      ,j=isp,isp+nnpts(i)-1)
    328:         write(13,1010) (1.0                 ,j=isp,isp+nnpts(i)-1)
    329:         write(13,1010) (q(2,isnode(j))      ,j=isp,isp+nnpts(i)-1)
    330:         write(13,1010) (q(3,isnode(j))      ,j=isp,isp+nnpts(i)-1)
    331:         write(13,1010) (q(4,isnode(j))      ,j=isp,isp+nnpts(i)-1)
    332:         write(13,1010) (q(1,isnode(j))/pinf ,j=isp,isp+nnpts(i)-1)

    334:         do 30 j=ist,ist+nntet(i)-1
    335:           i1 = f2ntn(j,1) - isp + 1
    336:           i2 = f2ntn(j,2) - isp + 1
    337:           i3 = f2ntn(j,3) - isp + 1

    339:           write(13,1020) i1,i2,i3,i3
    340:    30   continue

    342:         isp = isp + nnpts(i)
    343:         ist = ist + nntet(i)

    345:    10 continue
    346: !
    347: !   + start with viscous boundary surfaces
    348: !
    349:       isp   = 1
    350:       ist   = 1
    351:       do 70 i=1,nvbound
    352:         write(c4,"(i4)") i
    353:         if (i .ge.    0 .and. i .le.    9) ic = 4
    354:         if (i .ge.   10 .and. i .le.   99) ic = 3
    355:         if (i .ge.  100 .and. i .le.  999) ic = 2
    356:         if (i .ge. 1000 .and. i .le. 9999) ic = 1

    358:         write(13,1000) "nv." // c4(ic:4),nvpts(i),nvtet(i)

    360:         write(13,1010) (x(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    361:         write(13,1010) (y(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    362:         write(13,1010) (z(ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    363:         write(13,1010) (1.0                 ,j=isp,isp+nvpts(i)-1)
    364:         write(13,1010) (q(2,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    365:         write(13,1010) (q(3,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    366:         write(13,1010) (q(4,ivnode(j))      ,j=isp,isp+nvpts(i)-1)
    367:         write(13,1010) (q(1,ivnode(j))/pinf ,j=isp,isp+nvpts(i)-1)

    369:         do 90 j=ist,ist+nvtet(i)-1
    370:            i1 = f2ntv(j,1) - isp + 1
    371:            i2 = f2ntv(j,2) - isp + 1
    372:            i3 = f2ntv(j,3) - isp + 1

    374:            write(13,1020) i1,i2,i3,i3
    375:    90   continue

    377:         isp = isp + nvpts(i)
    378:         ist = ist + nvtet(i)

    380:    70 continue
    381: !
    382: !   + do the far-field boundary surfaces
    383: !
    384:       isp   = 1
    385:       ist   = 1
    386:       do 40 i=1,nfbound
    387:          write(c4,"(i4)") i
    388:          if (i .ge.    0 .and. i .le.    9) ic = 4
    389:          if (i .ge.   10 .and. i .le.   99) ic = 3
    390:          if (i .ge.  100 .and. i .le.  999) ic = 2
    391:          if (i .ge. 1000 .and. i .le. 9999) ic = 1

    393:          write(13,1000) "ff." // c4(ic:4),nfpts(i),nftet(i)

    395:          write(13,1010) (x(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
    396:          write(13,1010) (y(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
    397:          write(13,1010) (z(ifnode(j))        ,j=isp,isp+nfpts(i)-1)
    398:          write(13,1010) (1.0                 ,j=isp,isp+nfpts(i)-1)
    399:          write(13,1010) (q(2,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
    400:          write(13,1010) (q(3,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
    401:          write(13,1010) (q(4,ifnode(j))      ,j=isp,isp+nfpts(i)-1)
    402:          write(13,1010) (q(1,ifnode(j))/pinf ,j=isp,isp+nfpts(i)-1)

    404:          do 60 j=ist,ist+nftet(i)-1
    405:             i1 = f2ntf(j,1) - isp + 1
    406:             i2 = f2ntf(j,2) - isp + 1
    407:             i3 = f2ntf(j,3) - isp + 1

    409:             write(13,1020) i1,i2,i3,i3
    410:    60    continue

    412:          isp = isp + nfpts(i)
    413:          ist = ist + nftet(i)

    415:    40  continue

    417: !
    418: !      End of subroutine TECFLO
    419: !
    420:  1000  format('ZONE T="',a,'", Z=0, I=',i6,', J=',i6,', F=FEBLOCK')
    421:  1010  format(1P10E13.5)
    422:  1020  format(4I10)

    424:        return
    425:        end

    427: !================================ FORCE  =============================72
    428: !
    429: !  Calculates the forces
    430: !  Modified - D. K. Kaushik (1/15/97)
    431: !  Added new parameters - clift, cdrag, cmom, irank, nvertices
    432: !
    433: !=====================================================================72
    434:       subroutine FORCE(nnodes,nedge,                                    &
    435:      &                 isnode,ivnode,                                   &
    436:      &                 nnfacet,f2ntn,nnbound,                           &
    437:      &                 nvfacet,f2ntv,nvbound,                           &
    438:      &                 evec,qvec,xyz,sface_bit,vface_bit,               &
    439:      &                 clift, cdrag, cmom,irank,nvertices)
    440:       implicit none
    441: #include <petsc/finclude/petscsys.h>
    442:       PetscScalar  title(20),beta,alpha
    443:       PetscScalar  Re,dt,tot,res0,resc
    444:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    445:       PetscScalar gamma,gm1,gp1,gm1g,gp1g,ggm1
    446:       PetscScalar p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    447:       PetscScalar rms(1001),clw(1001),cdw(1001),                             &
    448:      &       cmw(1001),xres(1001)
    449:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                 &
    450:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    451:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
    452:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    453:       common/history/rms,clw,cdw,cmw,xres

    455:       integer nnodes,nedge,nnfacet,nvfacet,                             &
    456:      &        nnbound,nvbound,                                          &
    457:      &        irank,nvertices
    458:       integer isnode(1),ivnode(1)
    459:       integer f2ntn(nnfacet,4)
    460:       integer f2ntv(nvfacet,4)
    461:       integer sface_bit(nnfacet), vface_bit(nvfacet)

    463:       PetscScalar clift, cdrag, cmom
    464:       PetscScalar clloc,cdloc,cmloc
    465:       PetscScalar clglo,cdglo,cmglo
    466:       PetscScalar cl_ghost, cd_ghost, cm_ghost
    467:       PetscScalar coef_loc(3), coef_glo(3)
    468:       PetscScalar pi,conv,csa,sna,xr,yr,zr
    469:       PetscScalar x1,y1,z1,x2,y2,z2,x3,y3
    470:       PetscScalar z3,xmid,ymid,zmid
    471:       PetscScalar p1,p2,p3,u1,u2,u3,v1,v2,v3
    472:       PetscScalar w1,w2,w3,ax,ay,az,bx,by,bz
    473:       PetscScalar xnorm,ynorm,znorm
    474:       PetscScalar dcx,dcy,dcz,cp,press
    475:       integer icount,n,node1,node2,node3,ierr

    477: #if defined(INTERLACING)
    478:       PetscScalar qvec(4,nvertices)
    479:       PetscScalar xyz(3,nvertices)
    480:       integer evec(2,nedge)
    481: #define qnod(i,j) qvec(i,j)
    482: #define x(i) xyz(1,i)
    483: #define y(i) xyz(2,i)
    484: #define z(i) xyz(3,i)
    485: #define eptr(j,i) evec(i,j)
    486: #else
    487:       PetscScalar qvec(nvertices,4)
    488:       PetscScalar xyz(nvertices,3)
    489:       integer evec(nedge,2)
    490: #define qnod(i,j) qvec(j,i)
    491: #define x(i) xyz(i,1)
    492: #define y(i) xyz(i,2)
    493: #define z(i) xyz(i,3)
    494: #define eptr(i,j) evec(i,j)
    495: #endif
    496: !
    497:       pi = 4.d0*atan(1.d0)
    498:       conv = 180.d0/pi
    499:       csa=cos(alpha/conv)
    500:       sna=sin(alpha/conv)
    501: !
    502: !  initialize forces to zero
    503: !
    504:       clw(ntt) = 0.0d0
    505:       cdw(ntt) = 0.0d0
    506:       cmw(ntt) = 0.0d0
    507:       cl_ghost = 0.0d0
    508:       cd_ghost = 0.0d0
    509:       cm_ghost = 0.0d0
    510:       clloc = 0.0d0
    511:       cdloc = 0.0d0
    512:       cmloc = 0.0d0
    513:       clglo = 0.0d0
    514:       cdglo = 0.0d0
    515:       cmglo = 0.0d0
    516:       coef_glo(1) = 0.0d0
    517:       coef_glo(2) = 0.0d0
    518:       coef_glo(3) = 0.0d0

    520:       xr = 0.25d0
    521:       yr = 0.0d0
    522:       zr = 0.0d0

    524:       do 30 n = 1, nnfacet
    525:                node1 = isnode(f2ntn(n,1))
    526:                node2 = isnode(f2ntn(n,2))
    527:                node3 = isnode(f2ntn(n,3))

    529:                x1    = x(node1)
    530:                y1    = y(node1)
    531:                z1    = z(node1)

    533:                x2    = x(node2)
    534:                y2    = y(node2)
    535:                z2    = z(node2)

    537:                x3 = x(node3)
    538:                y3 = y(node3)
    539:                z3 = z(node3)

    541:                ax = x2 - x1
    542:                ay = y2 - y1
    543:                az = z2 - z1

    545:                bx = x3 - x1
    546:                by = y3 - y1
    547:                bz = z3 - z1
    548: !
    549: !  norm points outward, away from grid interior.
    550: !  norm magnitude is area of surface triangle.
    551: !
    552:                xnorm =-0.5d0*(ay*bz - az*by)
    553:                ynorm = 0.5d0*(ax*bz - az*bx)
    554:                znorm =-0.5d0*(ax*by - ay*bx)

    556:                p1 = qnod(1,node1)
    557:                u1 = qnod(2,node1)
    558:                v1 = qnod(3,node1)
    559:                w1 = qnod(4,node1)

    561:                p2 = qnod(1,node2)
    562:                u2 = qnod(2,node2)
    563:                v2 = qnod(3,node2)
    564:                w2 = qnod(4,node2)

    566:                p3 = qnod(1,node3)
    567:                u3 = qnod(2,node3)
    568:                v3 = qnod(3,node3)
    569:                w3 = qnod(4,node3)

    571:                press = (p1 + p2 + p3)/3.0d0
    572:                cp    = 2.d0*(press - 1.0d0)
    573: !
    574:                dcx = cp*xnorm
    575:                dcy = cp*ynorm
    576:                dcz = cp*znorm

    578:                xmid = (x1 + x2 + x3)
    579:                ymid = (y1 + y2 + y3)
    580:                zmid = (z1 + z2 + z3)

    582:                if (sface_bit(n) .eq. 1) then
    583: #if defined(CFL3D_AXIS)

    585:                 clloc = clloc - dcx*sna     + dcz*csa
    586:                 cdloc = cdloc + dcx*csa + dcz*sna
    587:                 cmloc = cmloc                                           &
    588:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx

    590: #else
    591:                 clloc = clloc - dcx*sna + dcy*csa
    592:                 cdloc = cdloc + dcx*csa + dcy*sna
    593:                 cmloc = cmloc                                           &
    594:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
    595: #endif
    596:                endif

    598:  30    continue
    599: !
    600: ! Viscous boundary
    601: !
    602:        do 60 n = 1, nvfacet
    603:                node1 = ivnode(f2ntv(n,1))
    604:                node2 = ivnode(f2ntv(n,2))
    605:                node3 = ivnode(f2ntv(n,3))

    607:                x1    = x(node1)
    608:                y1    = y(node1)
    609:                z1    = z(node1)

    611:                x2    = x(node2)
    612:                y2    = y(node2)
    613:                z2    = z(node2)

    615:                x3 = x(node3)
    616:                y3 = y(node3)
    617:                z3 = z(node3)

    619:                ax = x2 - x1
    620:                ay = y2 - y1
    621:                az = z2 - z1

    623:                bx = x3 - x1
    624:                by = y3 - y1
    625:                bz = z3 - z1
    626: !
    627: !  norm points outward, away from grid interior.
    628: !  norm magnitude is area of surface triangle.
    629: !
    630:                xnorm =-0.5d0*(ay*bz - az*by)
    631:                ynorm = 0.5d0*(ax*bz - az*bx)
    632:                znorm =-0.5d0*(ax*by - ay*bx)

    634:                p1    = qnod(1,node1)
    635:                u1    = qnod(2,node1)
    636:                v1    = qnod(3,node1)
    637:                w1    = qnod(4,node1)

    639:                p2    = qnod(1,node2)
    640:                u2    = qnod(2,node2)
    641:                v2    = qnod(3,node2)
    642:                w2    = qnod(4,node2)

    644:                p3    = qnod(1,node3)
    645:                u3    = qnod(2,node3)
    646:                v3    = qnod(3,node3)
    647:                w3    = qnod(4,node3)

    649:                press = (p1 + p2 + p3)/3.0d0
    650:                cp    = 2.d0*(press -1.0d0)

    652:                dcx = cp*xnorm
    653:                dcy = cp*ynorm
    654:                dcz = cp*znorm

    656:                xmid = (x1 + x2 + x3)
    657:                ymid = (y1 + y2 + y3)
    658:                zmid = (z1 + z2 + z3)

    660:                if (vface_bit(n) .eq. 1) then
    661:                 clloc = clloc - dcx*sna + dcy*csa
    662:                 cdloc = cdloc + dcx*csa + dcy*sna
    663:                 cmloc = cmloc                                            &
    664:      &                     + (xmid - xr)*dcy - (ymid - yr)*dcx
    665:                endif

    667:  60   continue
    668: !
    669:       icount = 3
    670:       coef_loc(1) = clloc
    671:       coef_loc(2) = cdloc
    672:       coef_loc(3) = cmloc
    673:       call MPI_ALLREDUCE(coef_loc,coef_glo,icount,                       &
    674:      &                   MPIU_SCALAR, MPIU_SUM,                           &
    675:      &                   MPI_COMM_WORLD,ierr)
    676: !     call MPI_ALLREDUCE(cdloc,cdglo,icount,
    677: !    &                   MPIU_REAL_PRECISION, MPI_SUM,
    678: !    &                   MPI_COMM_WORLD,ierr)
    679: !     call MPI_ALLREDUCE(cmloc,cmglo,icount,
    680: !    &                   MPIU_REAL_PRECISION, MPI_SUM,
    681: !    &                   MPI_COMM_WORLD,ierr)
    682: !
    683:        clw(ntt) = coef_glo(1)
    684:        cdw(ntt) = coef_glo(2)
    685:        cmw(ntt) = coef_glo(3)
    686:        clift = clw(ntt)
    687:        cdrag = cdw(ntt)
    688:        cmom = cmw(ntt)
    689: !
    690: ! Update the timing information - Added by D. K. Kaushik (1/17/97)
    691:       xres(ntt) = tot
    692: !
    693: ! End of subroutine FORCE
    694: !
    695:       return
    696:       end

    699: !================================ DELTAT2 ============================72
    700: !
    701: ! Calculate a time step for each cell
    702: ! Note that this routine assumes conservative variables
    703: !
    704: !=====================================================================72
    705:       subroutine DELTAT2(nnodes,nedge,qvec,cdt,                          &
    706:      &                  xyz,vol,xyzn,evec,                               &
    707:      &                  sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,             &
    708:      &                  nsnode,nvnode,nfnode,isnode,ivnode,ifnode,       &
    709:      &                  LocalTS,irank,nvertices)
    710: !
    711:       implicit none
    712: #include <petsc/finclude/petscsys.h>
    713:       PetscScalar  title(20),beta,alpha
    714:       PetscScalar  Re,dt,tot,res0,resc
    715:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    716:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
    717:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                  &
    718:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    719:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
    720:       integer nnodes,nedge,nsnode,nvnode,nfnode,LocalTS,irank,nvertices
    721:       PetscScalar  vol(1)
    722:       PetscScalar  sxn(1),syn(1),szn(1)
    723:       PetscScalar  vxn(1),vyn(1),vzn(1)
    724:       PetscScalar  fxn(1),fyn(1),fzn(1)
    725:       PetscScalar  cdt(1)
    726:       integer isnode(1),ivnode(1),ifnode(1)
    727:       integer n,i,node1,node2,inode
    728:       integer ierr
    729:       PetscLogDouble flops
    730:       PetscScalar xnorm,ynorm,znorm,area
    731:       PetscScalar u1,v1,w1,u2,v2,w2,u,v,c,w
    732:       PetscScalar ubar,Vn,term
    733: #if defined(INTERLACING)
    734:       PetscScalar qvec(4,nvertices)
    735:       PetscScalar xyz(3,nvertices)
    736:       PetscScalar xyzn(4,nedge)
    737:       integer evec(2,nedge)
    738: #define qnod(i,j) qvec(i,j)
    739: #define x(i) xyz(1,i)
    740: #define y(i) xyz(2,i)
    741: #define z(i) xyz(3,i)
    742: #define xn(i) xyzn(1,i)
    743: #define yn(i) xyzn(2,i)
    744: #define zn(i) xyzn(3,i)
    745: #define rl(i) xyzn(4,i)
    746: #define eptr(j,i) evec(i,j)
    747: #else
    748:       PetscScalar qvec(nvertices,4)
    749:       PetscScalar xyz(nvertices,3)
    750:       PetscScalar xyzn(nedge,4)
    751:       integer evec(nedge,2)
    752: #define qnod(i,j) qvec(j,i)
    753: #define x(i) xyz(i,1)
    754: #define y(i) xyz(i,2)
    755: #define z(i) xyz(i,3)
    756: #define xn(i) xyzn(i,1)
    757: #define yn(i) xyzn(i,2)
    758: #define zn(i) xyzn(i,3)
    759: #define rl(i) xyzn(i,4)
    760: #define eptr(i,j) evec(i,j)
    761: #endif
    762: !
    763: !  If local time steping, loop over faces
    764: !  and calculate time step as cdt = V/(sum(|u.n| +c.area)
    765: !  This is time step for cfl=1. We will multiply by cfl number later
    766: !
    767: ! First loop over nodes and zero out cdt
    768: !
    769:       flops = 0.0
    770:       if (LocalTS .gt. 0) then
    771:          do 1000 i = 1,nvertices
    772:             cdt(i) = 0.0d0
    773:  1000    continue
    774: !
    775:          do 1020 n = 1, nedge
    776:                node1 = eptr(n,1)
    777:                node2 = eptr(n,2)
    778: !
    779: ! Get normal to face
    780: !
    781:                xnorm = xn(n)
    782:                ynorm = yn(n)
    783:                znorm = zn(n)
    784:                area  = rl(n)
    785: !
    786:                xnorm = xnorm*area
    787:                ynorm = ynorm*area
    788:                znorm = znorm*area
    789: !
    790: !/*
    791: ! xnorm = xnormal * area of face
    792: ! ynorm = ynormal * area of face
    793: ! znorm = znormal * area of face
    794: !*/
    795:                u1   = qnod(2,node1)
    796:                v1   = qnod(3,node1)
    797:                w1   = qnod(4,node1)

    799:                u2   = qnod(2,node2)
    800:                v2   = qnod(3,node2)
    801:                w2   = qnod(4,node2)
    802: !
    803: ! Get average values on face
    804: !
    805:                u    = 0.5d0*(u1 + u2)
    806:                v    = 0.5d0*(v1 + v2)
    807:                w    = 0.5d0*(w1 + w2)
    808:                ubar = xn(n)*u + yn(n)*v + zn(n)*w
    809:                c    = sqrt(ubar*ubar + beta)

    811:                term = abs(u*xnorm + v*ynorm + w*znorm) + c*area
    812:                cdt(node1) = cdt(node1) + term
    813:                cdt(node2) = cdt(node2) + term
    814: !
    815:  1020    continue
    816:          flops = flops + 27.0*nedge
    817: !
    818: ! Now loop over boundaries and close the contours
    819: !
    820:          do 1030 i = 1,nsnode
    821:             inode = isnode(i)
    822: !
    823: ! Get the normal
    824: !
    825:             xnorm = sxn(i)
    826:             ynorm = syn(i)
    827:             znorm = szn(i)
    828:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)

    830:             u   = qnod(2,inode)
    831:             v   = qnod(3,inode)
    832:             w   = qnod(4,inode)
    833:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
    834:             c   = sqrt(ubar*ubar + beta)
    835: !
    836:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
    837:             cdt(inode) = cdt(inode) + Vn
    838: !
    839:  1030    continue
    840:          flops = flops + 24.0*nsnode

    842: !
    843: ! Now viscous faces
    844: !
    845: !
    846:          do 1040 i = 1,nvnode
    847:             inode = ivnode(i)
    848: !
    849: ! Get the normal
    850: !
    851:             xnorm = vxn(i)
    852:             ynorm = vyn(i)
    853:             znorm = vzn(i)
    854:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
    855: !
    856: ! Because velocities are zero, c is just sqrt(beta)
    857: !
    858:             c  = sqrt(beta)
    859:             Vn = c*area

    861:             cdt(inode) = cdt(inode) + Vn
    862: !
    863:  1040    continue
    864:          flops = flops + 9.0*nvnode
    865: !
    866: ! Now far field
    867: !
    868:          do 1050 i = 1,nfnode
    869:             inode = ifnode(i)
    870: !
    871: ! Get the normal
    872: !
    873:             xnorm = fxn(i)
    874:             ynorm = fyn(i)
    875:             znorm = fzn(i)
    876:             area  = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)

    878:             u   = qnod(2,inode)
    879:             v   = qnod(3,inode)
    880:             w   = qnod(4,inode)
    881:             ubar= (u*xnorm + v*ynorm + w*znorm)/area
    882:             c = sqrt(ubar*ubar + beta)

    884:             Vn = abs(xnorm*u + ynorm*v + znorm*w) + c*area
    885:             cdt(inode) = cdt(inode) + Vn
    886:  1050    continue
    887:          flops = flops + 24.0*nfnode

    889:          do 1060 n = 1,nvertices
    890:             cdt(n) = vol(n)/cdt(n)
    891:  1060    continue
    892:          flops = flops + nvertices
    893:       else
    894: !
    895: ! If not doing local time stepping just set cdt=1
    896: !
    897:          do 1070 n = 1,nvertices
    898:             cdt(n) = 1.0d0
    899:  1070    continue
    900:       end if
    901:       call PetscLogFlops(flops,ierr)
    902: !
    903: ! End of subroutine DELTAT2
    904: !
    905:       return
    906:       end

    909: !================================= FLUX ======================================
    910: !
    911: ! Calculates the fluxes on the face and performs the flux balance
    912: !
    913: !=============================================================================
    914: #if defined(_OPENMP)
    915: #if defined(HAVE_EDGE_COLORING)
    916:       subroutine FLUX(nnodes, ncell, nedge,                                     &
    917:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
    918:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    919:      &                nnfacet,f2ntn,nnbound,                                    &
    920:      &                nvfacet,f2ntv,nvbound,                                    &
    921:      &                nffacet,f2ntf,nfbound,                                    &
    922:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
    923:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    924:      &                max_threads,                                              &
    925:      &                ncolor,ncount,                                            &
    926:      &                irank,nvertices)
    927: #elif defined(HAVE_REDUNDANT_WORK)
    928:       subroutine FLUX(nnodes, ncell, nedge,                                     &
    929:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
    930:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    931:      &                nnfacet,f2ntn,nnbound,                                    &
    932:      &                nvfacet,f2ntv,nvbound,                                    &
    933:      &                nffacet,f2ntf,nfbound,                                    &
    934:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
    935:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    936:      &                max_threads,                                              &
    937:      &                resd,                                                     &
    938:      &                irank,nvertices)
    939: #else
    940:       subroutine FLUX(nnodes, ncell, nedge,                                     &
    941:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
    942:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    943:      &                nnfacet,f2ntn,nnbound,                                    &
    944:      &                nvfacet,f2ntv,nvbound,                                    &
    945:      &                nffacet,f2ntf,nfbound,                                    &
    946:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
    947:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    948:      &                max_threads,                                              &
    949:      &                nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr,         &
    950:      &                irank,nvertices)
    951: #endif
    952: #else
    953:       subroutine FLUX(nnodes, ncell, nedge,                                     &
    954:      &                nsface, nvface, nfface, isface, ivface, ifface,           &
    955:      &                nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    956:      &                nnfacet,f2ntn,nnbound,                                    &
    957:      &                nvfacet,f2ntv,nvbound,                                    &
    958:      &                nffacet,f2ntf,nfbound,                                    &
    959:      &                grad,evec,qvec,xyz,resvec,xyzn,                           &
    960:      &                sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    961:      &                irank,nvertices)
    962: #endif
    963:       implicit none
    964: #include <petsc/finclude/petscsys.h>
    965:       PetscScalar  title(20),beta,alpha
    966:       PetscScalar  Re,dt,tot,res0,resc
    967:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    968:       PetscScalar  gamma,gm1,gp1,gm1g,gp1g,ggm1
    969:       PetscScalar  p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    970:       PetscScalar  rms(1001),clw(1001),cdw(1001)
    971:       PetscScalar  cmw(1001),xres(1001)
    972:       PetscScalar  cfl1,cfl2
    973:       integer nsmoth,iflim,itran,nbtran,jupdate,                                &
    974:      &        nstage,ncyct,iramp,nitfo
    975:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                         &
    976:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    977:       common/fluid/gamma,gm1,gp1,gm1g,gp1g,ggm1
    978:       common/ivals/p0,rho0,c0,u0,v0,w0,et0,h0,pt0
    979:       common/history/rms,clw,cdw,cmw,xres
    980:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,                 &
    981:      &             nstage,ncyct,iramp,nitfo
    982:       integer nnodes, ncell, nedge,                                             &
    983:      &        nsface, nvface, nfface,                                           &
    984:      &        nsnode, nvnode, nfnode,                                           &
    985:      &        nnfacet,nvfacet,nffacet,                                          &
    986:      &        nnbound,nvbound,nfbound,                                          &
    987:      &        irank,nvertices
    988:       PetscScalar  sxn(1),syn(1),szn(1)
    989:       PetscScalar  vxn(1),vyn(1),vzn(1)
    990:       PetscScalar  fxn(1),fyn(1),fzn(1)
    991:       PetscScalar  phi(nvertices,4)
    992:       integer isface(1),ivface(1),ifface(1)
    993:       integer isnode(1),ivnode(1),ifnode(1)
    994:       integer f2ntn(nnfacet,4)
    995:       integer f2ntv(nvfacet,4)
    996:       integer f2ntf(nffacet,4)
    997:       integer i,n,node1,node2,node3,inode
    998:       integer ierr,first_edge,last_edge
    999:       PetscLogDouble flops
    1000:       PetscScalar  xmean,ymean,zmean,xnorm,ynorm
    1001:       PetscScalar  znorm,rlen,area,dot,size
    1002:       PetscScalar  X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3
    1003:       PetscScalar  p1,p2,p3
    1004:       PetscScalar  rx,ry,rz,pL,uL,wL,vL
    1005:       PetscScalar  pR,uR,vR,wR
    1006:       PetscScalar  ubarL,ubarR,ubar,p,u,v,w,c
    1007:       PetscScalar  c2,c20
    1008:       PetscScalar  ubar0,pi,ui,vi,wi,unorm
    1009:       PetscScalar  phi1,phi2,phi3,phi4,phi5
    1010:       PetscScalar  phi6,phi7,phi8,phi9
    1011:       PetscScalar  eig1,eig2,eig3,eig4
    1012:       PetscScalar  dp,du,dv,dw
    1013:       PetscScalar  ti11,ti21,ti31,ti41
    1014:       PetscScalar  dv1,dv2,dv3,dv4
    1015:       PetscScalar  r11,r21,r31,r41
    1016:       PetscScalar  r12,r22,r32,r42
    1017:       PetscScalar  r13,r23,r33,r43
    1018:       PetscScalar  r14,r24,r34,r44
    1019:       PetscScalar  t1,t2,t3,t4
    1020:       PetscScalar  t11,t21,t31,t41
    1021:       PetscScalar  t12,t22,t32,t42
    1022:       PetscScalar  t13,t23,t33,t43
    1023:       PetscScalar  t14,t24,t34,t44
    1024:       PetscScalar  fluxp1,fluxp2,fluxp3,fluxp4
    1025:       PetscScalar  fluxm1,fluxm2,fluxm3,fluxm4
    1026:       PetscScalar  res1,res2,res3,res4
    1027:       PetscScalar  rhs1,rhs2,rhs3,rhs4
    1028:       PetscScalar  c68,c18,second
    1029:       PetscScalar  ax,ay,az,bx,by,bz
    1030:       PetscScalar  pa,pb,pc,ub,vb,wb
    1031: #if defined(INTERLACING)
    1032:       PetscScalar qvec(4,nvertices)
    1033:       PetscScalar grad(3,4,nvertices)
    1034:       PetscScalar resvec(4,nnodes)
    1035:       PetscScalar xyz(3,nvertices)
    1036:       PetscScalar xyzn(4,nedge)
    1037:       integer evec(2,nedge)
    1038: #define qnod(i,j) qvec(i,j)
    1039: #define res(i,j) resvec(i,j)
    1040: #define gradx(i,j) grad(1,i,j)
    1041: #define grady(i,j) grad(2,i,j)
    1042: #define gradz(i,j) grad(3,i,j)
    1043: #define x(i) xyz(1,i)
    1044: #define y(i) xyz(2,i)
    1045: #define z(i) xyz(3,i)
    1046: #define xn(i) xyzn(1,i)
    1047: #define yn(i) xyzn(2,i)
    1048: #define zn(i) xyzn(3,i)
    1049: #define ra(i) xyzn(4,i)
    1050: #define eptr(j,i) evec(i,j)
    1051: #else
    1052:       PetscScalar qvec(nvertices,4)
    1053:       PetscScalar grad(nvertices,4,3)
    1054:       PetscScalar resvec(nnodes,4)
    1055:       PetscScalar xyz(nvertices,3)
    1056:       PetscScalar xyzn(nedge,4)
    1057:       integer evec(nedge,2)
    1058: #define qnod(i,j) qvec(j,i)
    1059: #define res(i,j) resvec(j,i)
    1060: #define gradx(i,j) grad(j,i,1)
    1061: #define grady(i,j) grad(j,i,2)
    1062: #define gradz(i,j) grad(j,i,3)
    1063: #define x(i) xyz(i,1)
    1064: #define y(i) xyz(i,2)
    1065: #define z(i) xyz(i,3)
    1066: #define xn(i) xyzn(i,1)
    1067: #define yn(i) xyzn(i,2)
    1068: #define zn(i) xyzn(i,3)
    1069: #define ra(i) xyzn(i,4)
    1070: #define eptr(i,j) evec(i,j)
    1071: #endif
    1072: #if defined(_OPENMP)
    1073:        integer my_thread_id,omp_get_thread_num,                             &
    1074:      &         max_threads,omp_get_num_threads,                             &
    1075:      &         first_node,last_node,chunk
    1076: #if defined(HAVE_EDGE_COLORING)
    1077:        integer ncolor,ncount(1)
    1078: #elif defined(HAVE_REDUNDANT_WORK)
    1079:        PetscScalar resd(4,nnodes)
    1080: #undef res
    1081: #define res(a,b) resd(a,b)
    1082: #else
    1083:        integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads),         &
    1084:      &     edge_thr(2,nedgeAllThr)
    1085:        PetscScalar xyzn_thr(4,nedgeAllThr)
    1086: #undef eptr
    1087: #define eptr(a,b) edge_thr(b,a)
    1088: #undef xn
    1089: #undef yn
    1090: #undef zn
    1091: #undef ra
    1092: #define xn(i) xyzn_thr(1,i)
    1093: #define yn(i) xyzn_thr(2,i)
    1094: #define zn(i) xyzn_thr(3,i)
    1095: #define ra(i) xyzn_thr(4,i)
    1096: #endif
    1097: #endif
    1098: !      logging variables
    1099:        integer flux_event, fun3d_class, flag
    1100:        character * 16 flux_label
    1101:        data flag/-1/,flux_label/'FluxEvaluation  '/
    1102:        save flux_event, fun3d_class, flag, flux_label
    1103:        if (flag .eq. -1) then
    1104:           call PetscClassIdRegister('PetscFun3d',fun3d_class,ierr)
    1105:           call PetscLogEventRegister(flux_label,                        &
    1106:      &         fun3d_class,flux_event,ierr)
    1107:           flag = 1
    1108:        endif
    1109:        call PetscLogEventBegin(flux_event,ierr)
    1110:        flops = 0.0
    1111:        second = 1.0d0
    1112: !      if (ntt.le.nitfo)second = 0.0
    1113: !      print *, 'Second is' , second
    1114: #if defined(_OPENMP)
    1115: #if defined(HAVE_EDGE_COLORING)
    1116:        first_edge = 1
    1117:        do i = 1,ncolor
    1118:              last_edge = first_edge + ncount(i) - 1
    1119: #elif defined(HAVE_REDUNDANT_WORK)
    1120:        do i = 1,nnodes
    1121:           do j = 1,4
    1122:              resd(j,i) = 0.0
    1123:           enddo
    1124:        enddo
    1125:        first_edge = 1
    1126:        last_edge = nedge
    1127: #endif
    1128: #else
    1129:        first_edge = 1
    1130:        last_edge = nedge
    1131: #endif
    1132: !
    1133: ! Loop over all the faces and calculate flux
    1134: !
    1135: !$omp  parallel default(shared) private(n,node1,node2,xmean,ymean,           &
    1136: !$omp& zmean,xnorm,ynorm,znorm,rlen,dot,x1,y1,z1,size,x2,y2,z2,rx,           &
    1137: !$omp& ry,rz,pL,uL,vL,wL,ubarL,pR,uR,vR,wR,ubarR,p,u,v,w,ubar,phi1,          &
    1138: !$omp& phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,c2,c,eig1,eig2,eig3,          &
    1139: !$omp& eig4,dp,du,dv,dw,ti11,ti21,ti31,ti41,dv1,dv2,dv3,dv4,r11,r21,         &
    1140: !$omp& r31,r41,r12,r22,r32,r42,r13,r23,r33,r43,r14,r24,r34,r44,t1,t2,        &
    1141: !$omp& t3,t4,fluxp1,fluxp2,fluxp3,fluxp4,fluxm1,fluxm2,fluxm3,fluxm4,        &
    1142: !$omp& res1,res2,res3,res4,c68,c18,inode)                                    &
    1143: !$omp& private(x3,y3,z3,p1,p2,p3,ax,ay,az,bx,by,bz,pa,pb,pc,area,c20,        &
    1144: !$omp& c0,ubar0,t11,t21,t31,t41,t12,t22,t32,t42,t13,t23,t33,t43,t14,         &
    1145: !$omp& t24,t34,t44,pi,ui,vi,wi,unorm,rhs1,rhs2,rhs3,rhs4,ub,vb,wb,node3)     &
    1146: !$omp& private(my_thread_id,chunk)                                           &
    1147: !$omp& reduction(+:flops)                                                    &
    1148: #if defined(_OPENMP)
    1149: #if defined(HAVE_EDGE_COLORING)
    1150: #elif defined(HAVE_REDUNDANT_WORK)
    1151: !$omp& firstprivate(resd)
    1152: #else
    1153: !$omp& private(first_edge,last_edge)
    1154: #endif
    1155: #endif
    1156: !!$omp1 shared(evec,qvec,resvec,xyz,xyzn,grad,nedge,nnodes,second,flops,
    1157: !!$omp2 title,beta,alpha,Re,dt,tot,res0,resc,ntt,mseq,ivisc,irest,icyc,
    1158: !!$omp3 ihane,ntturb,gamma,gm1,gp1,gm1g,gp1g,ggm1,p0,rho0,c0,u0,v0,w0,
    1159: !!$omp4 et0,h0,pt0,rms,clw,cdw,cmw,xres,cfl1,cfl2,nsmoth,iflim,
    1160: !!$omp5 itran,nbtran,jupdate,nstage,ncyct,iramp,nitfo)
    1161: !!$      my_thread_id = omp_get_thread_num()
    1162: !!$omp   critical
    1163: !!$      print *, 'My thread id on processor ',irank,' is ',my_thread_id
    1164: !!$omp   end critical
    1165: #if defined(_OPENMP)
    1166: #if defined(HAVE_EDGE_COLORING)
    1167: !$omp do
    1168: #elif defined(HAVE_REDUNDANT_WORK)
    1169: !$omp do
    1170: #else
    1171: !     tot_threads = omp_get_num_threads()
    1172:       my_thread_id = omp_get_thread_num()
    1173: !     chunk = nnodes/tot_threads
    1174: !     first_node = my_thread_id*chunk+1
    1175: !     if (my_thread_id .eq. (size-1)) then
    1176: !        last_node = nnodes
    1177: !     else
    1178: !        last_node = (my_thread_id+1)*chunk
    1179: !     endif
    1180:       first_edge = nedge_thr(my_thread_id)
    1181:       last_edge = nedge_thr(my_thread_id+1)-1

    1183: !      print *, 'On thread ', my_thread_id,', first_edge = ',first_edge,        &
    1184: !     &         ', last_edge = ',last_edge

    1186: #endif
    1187: #else
    1188: #endif
    1189:       do n = first_edge, last_edge
    1190:         node1 = eptr(n,1)
    1191:         node2 = eptr(n,2)
    1192: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
    1193:         if ((part_thr(node1) .eq. my_thread_id)                                 &
    1194:      &       .or.(part_thr(node2) .eq. my_thread_id)) then
    1195: #endif
    1196: !
    1197: ! Calculate unit normal to face and length of face
    1198: !
    1199:           xmean = .5d0*(x(node1) + x(node2))
    1200:           ymean = .5d0*(y(node1) + y(node2))
    1201:           zmean = .5d0*(z(node1) + z(node2))
    1202:           xnorm  = xn(n)
    1203:           ynorm  = yn(n)
    1204:           znorm  = zn(n)
    1205:           rlen   = ra(n)
    1206: !
    1207: ! Now lets get our other 2 vectors
    1208: ! For first vector, use {1,0,0} and subtract off the component
    1209: ! in the direction of the face normal. If the inner product of
    1210: ! {1,0,0} is close to unity, use {0,1,0}
    1211: !
    1212:          dot = xnorm
    1213:          if (abs(dot).lt.0.95d0)then
    1214:           X1 = 1. - dot*xnorm
    1215:           Y1 =    - dot*ynorm
    1216:           Z1 =    - dot*znorm
    1217:          else
    1218:           dot = ynorm
    1219:           X1 =    - dot*xnorm
    1220:           Y1 = 1.0d0 - dot*ynorm
    1221:           Z1 =    - dot*znorm
    1222:          end if
    1223: !
    1224: ! Normalize the first vector
    1225: !
    1226:          size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
    1227:          X1 = X1/size
    1228:          Y1 = Y1/size
    1229:          Z1 = Z1/size
    1230: !
    1231: ! Take cross-product of normal and V1 to get V2
    1232: !
    1233:          X2 = ynorm*Z1 - znorm*Y1
    1234:          Y2 = znorm*X1 - xnorm*Z1
    1235:          Z2 = xnorm*Y1 - ynorm*X1
    1236: !
    1237: ! Get variables on "left" and "right" side of face
    1238: !
    1239:          rx     = second*(xmean - x(node1))
    1240:          ry     = second*(ymean - y(node1))
    1241:          rz     = second*(zmean - z(node1))
    1242:          pL = qnod(1,node1) + gradx(1,node1)*rx                           &
    1243:      &                       + grady(1,node1)*ry                           &
    1244:      &                       + gradz(1,node1)*rz
    1245:          uL = qnod(2,node1) + gradx(2,node1)*rx                           &
    1246:      &                       + grady(2,node1)*ry                           &
    1247:      &                       + gradz(2,node1)*rz
    1248:          vL = qnod(3,node1) + gradx(3,node1)*rx                           &
    1249:      &                       + grady(3,node1)*ry                           &
    1250:      &                       + gradz(3,node1)*rz
    1251:          wL = qnod(4,node1) + gradx(4,node1)*rx                           &
    1252:      &                       + grady(4,node1)*ry                           &
    1253:      &                       + gradz(4,node1)*rz

    1255:          ubarL  = xnorm*uL + ynorm*vL + znorm*wL
    1256: !        c2L = ubarl*ubarL + beta
    1257: !        cL  = sqrt(c2L)
    1258: !
    1259:          rx     = second*(xmean - x(node2))
    1260:          ry     = second*(ymean - y(node2))
    1261:          rz     = second*(zmean - z(node2))
    1262:          pR   = qnod(1,node2) + gradx(1,node2)*rx                         &
    1263:      &                         + grady(1,node2)*ry                         &
    1264:      &                         + gradz(1,node2)*rz
    1265:          uR   = qnod(2,node2) + gradx(2,node2)*rx                         &
    1266:      &                         + grady(2,node2)*ry                         &
    1267:      &                         + gradz(2,node2)*rz
    1268:          vR   = qnod(3,node2) + gradx(3,node2)*rx                         &
    1269:      &                         + grady(3,node2)*ry                         &
    1270:      &                         + gradz(3,node2)*rz
    1271:          wR   = qnod(4,node2) + gradx(4,node2)*rx                         &
    1272:      &                         + grady(4,node2)*ry                         &
    1273:      &                         + gradz(4,node2)*rz
    1274:          ubarR  = xnorm*uR + ynorm*vR + znorm*wR
    1275: !        c2R = ubarR*ubarR + beta
    1276: !        cR  = sqrt(c2R)
    1277: !
    1278: ! Compute averages
    1279: !
    1280:           p = .5d0*(pL + pR)
    1281:           u = .5d0*(uL + uR)
    1282:           v = .5d0*(vL + vR)
    1283:           w = .5d0*(wL + wR)

    1285:           ubar  = xnorm*u + ynorm*v + znorm*w
    1286:           phi1  = xnorm*beta + u*ubar
    1287:           phi2  = ynorm*beta + v*ubar
    1288:           phi3  = znorm*beta + w*ubar
    1289:           phi4  = Y2*phi3 - Z2*phi2
    1290:           phi5  = Z2*phi1 - X2*phi3
    1291:           phi6  = X2*phi2 - Y2*phi1
    1292:           phi7  = Z1*phi2 - Y1*phi3
    1293:           phi8  = X1*phi3 - Z1*phi1
    1294:           phi9  = Y1*phi1 - X1*phi2
    1295:           c2    = ubar*ubar + beta
    1296:           c     = sqrt(c2)
    1297: !
    1298: ! Now compute eigenvalues, eigenvectors, and strengths
    1299: !
    1300:           eig1 = abs(ubar)
    1301:           eig2 = abs(ubar)
    1302:           eig3 = abs(ubar + c)
    1303:           eig4 = abs(ubar - c)
    1304: !
    1305:           dp = pr - pl
    1306:           du = ur - ul
    1307:           dv = vr - vl
    1308:           dw = wr - wl
    1309: !
    1310: ! Components of T(inverse) (I will divide by c2 later)
    1311: !
    1312:           ti11 = -(u*phi4 + v*phi5 + w*phi6)/beta
    1313:           ti21 = -(u*phi7 + v*phi8 + w*phi9)/beta
    1314:           ti31 =  .5d0*(c-ubar)/beta
    1315:           ti41 = -.5d0*(c+ubar)/beta
    1316: !
    1317: ! jumps (T(inverse)*dq)
    1318:           dv1 = (ti11*dp + phi4*du + phi5*dv + phi6*dw)/c2
    1319:           dv2 = (ti21*dp + phi7*du + phi8*dv + phi9*dw)/c2
    1320:           dv3 = .5d0*(2.d0*ti31*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
    1321:           dv4 = .5d0*(2.d0*ti41*dp + xnorm*du + ynorm*dv + znorm*dw)/c2
    1322: !
    1323: ! Now get elements of T (call it r for now)
    1324: !
    1325:           r11 = 0.0d0
    1326:           r21 = X1
    1327:           r31 = Y1
    1328:           r41 = Z1

    1330:           r12 = 0.0d0
    1331:           r22 = X2
    1332:           r32 = Y2
    1333:           r42 = Z2

    1335:           r13 = c*beta
    1336:           r23 = xnorm*beta + u*(ubar + c)
    1337:           r33 = ynorm*beta + v*(ubar + c)
    1338:           r43 = znorm*beta + w*(ubar + c)
    1339: !
    1340:           r14 = -c*beta
    1341:           r24 = xnorm*beta + u*(ubar - c)
    1342:           r34 = ynorm*beta + v*(ubar - c)
    1343:           r44 = znorm*beta + w*(ubar - c)
    1344: !
    1345: ! Calculate T* |lambda| *T(inverse)
    1346: !
    1347:           t1 = eig1*r11*dv1 + eig2*r12*dv2 + eig3*r13*dv3 + eig4*r14*dv4
    1348:           t2 = eig1*r21*dv1 + eig2*r22*dv2 + eig3*r23*dv3 + eig4*r24*dv4
    1349:           t3 = eig1*r31*dv1 + eig2*r32*dv2 + eig3*r33*dv3 + eig4*r34*dv4
    1350:           t4 = eig1*r41*dv1 + eig2*r42*dv2 + eig3*r43*dv3 + eig4*r44*dv4
    1351: !
    1352: ! Modify to calculate .5(fl +fr) from nodes
    1353: ! instead of extrapolated ones
    1354: !
    1355: !         pL    = qnod(1,node1)
    1356: !         uL    = qnod(2,node1)
    1357: !         vL    = qnod(3,node1)
    1358: !         wL    = qnod(4,node1)
    1359: !         ubarL = xnorm*uL + ynorm*vL + znorm*wL
    1360: !
    1361:           fluxp1 = rlen*beta*ubarL
    1362:           fluxp2 = rlen*(uL*ubarL + xnorm*pL)
    1363:           fluxp3 = rlen*(vL*ubarL + ynorm*pL)
    1364:           fluxp4 = rlen*(wL*ubarL + znorm*pL)
    1365: !
    1366: ! Now the right side
    1367: !
    1368: !         pR    = qnod(1,node2)
    1369: !         uR    = qnod(2,node2)
    1370: !         vR    = qnod(3,node2)
    1371: !         wR    = qnod(4,node2)
    1372: !         ubarR = xnorm*uR + ynorm*vR + znorm*wR
    1373: !
    1374:           fluxm1 = rlen*beta*ubarR
    1375:           fluxm2 = rlen*(uR*ubarR + xnorm*pR)
    1376:           fluxm3 = rlen*(vR*ubarR + ynorm*pR)
    1377:           fluxm4 = rlen*(wR*ubarR + znorm*pR)
    1378: !
    1379:           res1 = 0.5d0*(fluxp1 + fluxm1 - rlen*t1)
    1380:           res2 = 0.5d0*(fluxp2 + fluxm2 - rlen*t2)
    1381:           res3 = 0.5d0*(fluxp3 + fluxm3 - rlen*t3)
    1382:           res4 = 0.5d0*(fluxp4 + fluxm4 - rlen*t4)
    1383: !
    1384:           flops = flops + 318.0

    1386: !
    1387: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
    1388:           if ((part_thr(node1).eq.my_thread_id)                                 &
    1389:      &        .and.(node1.le.nnodes)) then
    1390: #else
    1391:           if (node1.le.nnodes) then
    1392: #endif
    1393:            res(1,node1) = res(1,node1) + res1
    1394:            res(2,node1) = res(2,node1) + res2
    1395:            res(3,node1) = res(3,node1) + res3
    1396:            res(4,node1) = res(4,node1) + res4
    1397:            flops = flops + 4.0
    1398:           endif
    1399: !
    1400: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
    1401:           if ((part_thr(node2).eq.my_thread_id)                                 &
    1402:      &        .and.(node2.le.nnodes)) then
    1403: #else
    1404:           if (node2.le.nnodes) then
    1405: #endif
    1406:            res(1,node2) = res(1,node2) - res1
    1407:            res(2,node2) = res(2,node2) - res2
    1408:            res(3,node2) = res(3,node2) - res3
    1409:            res(4,node2) = res(4,node2) - res4
    1410:            flops = flops + 4.0
    1411:           endif
    1412: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK)&& !defined(HAVE_EDGE_COLORING)
    1413:          endif
    1414: #endif
    1415: !
    1416:       enddo
    1417: #if defined(_OPENMP)
    1418: #if defined(HAVE_EDGE_COLORING)
    1419: !$omp enddo
    1420: !$omp end parallel
    1421:        first_edge = first_edge + ncount(i)
    1422:       enddo
    1423: #elif defined(HAVE_REDUNDANT_WORK)
    1424: !$omp enddo
    1425: !$omp critical (residual_update)
    1426:        do n = 1,nnodes
    1427:           resvec(1,n)=resvec(1,n)+res(1,n)
    1428:           resvec(2,n)=resvec(2,n)+res(2,n)
    1429:           resvec(3,n)=resvec(3,n)+res(3,n)
    1430:           resvec(4,n)=resvec(4,n)+res(4,n)
    1431:        enddo
    1432:        flops = flops + 4.0*nnodes
    1433: !$omp end critical (residual_update)
    1434: !$omp end parallel
    1435: #undef res
    1436: #define res(a,b) resvec(a,b)
    1437: #else
    1438: !$omp end parallel
    1439: #endif
    1440: #else
    1441: #endif
    1442:       call PetscLogFlops(flops,ierr)
    1443:       call PetscLogEventEnd(flux_event,ierr)
    1444: !     c68 = 6./8.
    1445: !     c18 = 1./8.
    1446:       c68 = 0.75d0
    1447:       c18 = 0.125d0
    1448: !
    1449: ! Close contour over the boundaries
    1450: ! First do inviscid faces
    1451: !
    1452:       do n = 1, nnfacet
    1453:                node1 = isnode(f2ntn(n,1))
    1454:                node2 = isnode(f2ntn(n,2))
    1455:                node3 = isnode(f2ntn(n,3))

    1457:                x1 = x(node1)
    1458:                y1 = y(node1)
    1459:                z1 = z(node1)
    1460:                p1 = qnod(1,node1)

    1462:                x2 = x(node2)
    1463:                y2 = y(node2)
    1464:                z2 = z(node2)
    1465:                p2 = qnod(1,node2)

    1467:                x3 = x(node3)
    1468:                y3 = y(node3)
    1469:                z3 = z(node3)
    1470:                p3 = qnod(1,node3)

    1472:                ax = x2 - x1
    1473:                ay = y2 - y1
    1474:                az = z2 - z1

    1476:                bx = x3 - x1
    1477:                by = y3 - y1
    1478:                bz = z3 - z1
    1479: !
    1480: ! Normal points away from grid interior.
    1481: ! Magnitude is 1/3 area of surface triangle.
    1482: !
    1483:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
    1484:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
    1485:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

    1487:                pa = c68*p1 + c18*(p2 + p3)
    1488:                pb = c68*p2 + c18*(p3 + p1)
    1489:                pc = c68*p3 + c18*(p1 + p2)
    1490: !
    1491:                flops = flops + 35.0
    1492:           if (node1.le.nnodes) then
    1493:                res(2,node1) = res(2,node1) + xnorm*pa
    1494:                res(3,node1) = res(3,node1) + ynorm*pa
    1495:                res(4,node1) = res(4,node1) + znorm*pa
    1496:                flops = flops + 6.0
    1497:           endif

    1499:           if (node2.le.nnodes) then
    1500:                res(2,node2) = res(2,node2) + xnorm*pb
    1501:                res(3,node2) = res(3,node2) + ynorm*pb
    1502:                res(4,node2) = res(4,node2) + znorm*pb
    1503:                flops = flops + 6.0
    1504:           endif

    1506:           if (node3.le.nnodes) then
    1507:                res(2,node3) = res(2,node3) + xnorm*pc
    1508:                res(3,node3) = res(3,node3) + ynorm*pc
    1509:                res(4,node3) = res(4,node3) + znorm*pc
    1510:                flops = flops + 6.0
    1511:           endif

    1513:       enddo
    1514: !
    1515: ! Now viscous faces
    1516: !
    1517:       do n = 1, nvfacet
    1518:                node1 = ivnode(f2ntv(n,1))
    1519:                node2 = ivnode(f2ntv(n,2))
    1520:                node3 = ivnode(f2ntv(n,3))

    1522:                x1 = x(node1)
    1523:                y1 = y(node1)
    1524:                z1 = z(node1)
    1525:                p1 = qnod(1,node1)

    1527:                x2 = x(node2)
    1528:                y2 = y(node2)
    1529:                z2 = z(node2)
    1530:                p2 = qnod(1,node2)

    1532:                x3 = x(node3)
    1533:                y3 = y(node3)
    1534:                z3 = z(node3)
    1535:                p3 = qnod(1,node3)

    1537:                ax = x2 - x1
    1538:                ay = y2 - y1
    1539:                az = z2 - z1

    1541:                bx = x3 - x1
    1542:                by = y3 - y1
    1543:                bz = z3 - z1
    1544: !
    1545: ! norm point away from grid interior.
    1546: ! norm magnitude is 1/3 area of surface triangle.
    1547: !
    1548:                xnorm =-0.5d0*(ay*bz - az*by)/3.d0
    1549:                ynorm = 0.5d0*(ax*bz - az*bx)/3.d0
    1550:                znorm =-0.5d0*(ax*by - ay*bx)/3.d0

    1552:                pa = c68*p1 + c18*(p2 + p3)
    1553:                pb = c68*p2 + c18*(p3 + p1)
    1554:                pc = c68*p3 + c18*(p1 + p2)
    1555: !
    1556:                flops = flops + 35.0
    1557: !
    1558:           if (node1.le.nnodes) then
    1559:                 res(2,node1) = res(2,node1) + xnorm*pa
    1560:                 res(3,node1) = res(3,node1) + ynorm*pa
    1561:                 res(4,node1) = res(4,node1) + znorm*pa
    1562:                 flops = flops + 6.0

    1564:           endif

    1566:           if (node2.le.nnodes) then
    1567:                 res(2,node2) = res(2,node2) + xnorm*pb
    1568:                 res(3,node2) = res(3,node2) + ynorm*pb
    1569:                 res(4,node2) = res(4,node2) + znorm*pb
    1570:                 flops = flops + 6.0
    1571:           endif

    1573:           if (node3.le.nnodes) then
    1574:                 res(2,node3) = res(2,node3) + xnorm*pc
    1575:                 res(3,node3) = res(3,node3) + ynorm*pc
    1576:                 res(4,node3) = res(4,node3) + znorm*pc
    1577:                 flops = flops + 6.0
    1578:           endif

    1580:       enddo
    1581: !
    1582: ! The next section of code is for when you dont care about
    1583: ! preserving linear data on boundary. Also, doing this
    1584: ! matches the left hand side when not doing Newton-Krylov
    1585: ! Usually just go around unless you are just experimenting
    1586: !
    1587:       goto 1025
    1588: !
    1589: ! Loop over the boundaries
    1590: ! First do inviscid nodes
    1591: !
    1592: !     do 1020 i = 1,nsnode
    1593: !       inode = isnode(i)
    1594: !
    1595: !       xnorm = sxn(i)
    1596: !       ynorm = syn(i)
    1597: !       znorm = szn(i)
    1598: !
    1599: !       p = qnod(1,inode)
    1600: !
    1601: !       if (inode .le. nnodes) then
    1602: !        res(2,inode) = res(2,inode) + xnorm*p
    1603: !        res(3,inode) = res(3,inode) + ynorm*p
    1604: !        res(4,inode) = res(4,inode) + znorm*p
    1605: !       endif
    1606: !
    1607: !1020 continue
    1608: !
    1609: ! Now viscous nodes
    1610: !
    1611: !     do 1030 i = 1,nvnode
    1612: !       inode = ivnode(i)
    1613: !
    1614: !       xnorm = vxn(i)
    1615: !       ynorm = vyn(i)
    1616: !       znorm = vzn(i)
    1617: !
    1618: !       p = qnod(1,inode)
    1619: !
    1620: !       if (inode .le. nnodes) then
    1621: !        res(2,inode) = res(2,inode) + xnorm*p
    1622: !        res(3,inode) = res(3,inode) + ynorm*p
    1623: !        res(4,inode) = res(4,inode) + znorm*p
    1624: !       endif
    1625: !
    1626: !1030 continue
    1627: !
    1628:  1025 continue
    1629: !
    1630: ! Now do far-field
    1631: !
    1632:        do i = 1,nfnode
    1633:          inode   = ifnode(i)
    1634: !
    1635: !/*
    1636: ! Get normal and "other" 2 vectors. Remember that fxn,fyn and fzn
    1637: ! has the magnitude of the face contained in it.
    1638: !*/
    1639:          xnorm   = fxn(i)
    1640:          ynorm   = fyn(i)
    1641:          znorm   = fzn(i)
    1642:          area    = sqrt(xnorm*xnorm + ynorm*ynorm + znorm*znorm)
    1643:          xnorm   = xnorm/area
    1644:          ynorm   = ynorm/area
    1645:          znorm   = znorm/area
    1646: !
    1647: ! Now lets get our other 2 vectors
    1648: ! For first vector, use {1,0,0} and subtract off the component
    1649: ! in the direction of the face normal. If the inner product of
    1650: ! {1,0,0} is close to unity, use {0,1,0}
    1651: !
    1652:          dot = xnorm
    1653:          if (abs(dot).lt.0.95d0)then
    1654:           X1 = 1.d0 - dot*xnorm
    1655:           Y1 =    - dot*ynorm
    1656:           Z1 =    - dot*znorm
    1657:          else
    1658:           dot = ynorm
    1659:           X1 =    - dot*xnorm
    1660:           Y1 = 1.d0 - dot*ynorm
    1661:           Z1 =    - dot*znorm
    1662:          end if
    1663: !
    1664: ! Normalize the first vector (V1)
    1665: !
    1666:          size = sqrt(X1*X1 + Y1*Y1 + Z1*Z1)
    1667:          X1 = X1/size
    1668:          Y1 = Y1/size
    1669:          Z1 = Z1/size
    1670: !
    1671: ! Take cross-product of normal with V1 to get V2
    1672: !
    1673:          X2 = ynorm*Z1 - znorm*Y1
    1674:          Y2 = znorm*X1 - xnorm*Z1
    1675:          Z2 = xnorm*Y1 - ynorm*X1

    1677: !
    1678: ! Calculate elements of T and T(inverse) evaluated at freestream
    1679: !
    1680:          ubar0 = xnorm*u0 + ynorm*v0 + znorm*w0
    1681:          c20   = ubar0*ubar0 + beta
    1682:          c0    = sqrt(c20)
    1683:          phi1  = xnorm*beta + u0*ubar0
    1684:          phi2  = ynorm*beta + v0*ubar0
    1685:          phi3  = znorm*beta + w0*ubar0
    1686:          phi4  = Y2*phi3 - Z2*phi2
    1687:          phi5  = Z2*phi1 - X2*phi3
    1688:          phi6  = X2*phi2 - Y2*phi1
    1689:          phi7  = Z1*phi2 - Y1*phi3
    1690:          phi8  = X1*phi3 - Z1*phi1
    1691:          phi9  = Y1*phi1 - X1*phi2

    1693:          t11 = 0.0d0
    1694:          t21 = X1
    1695:          t31 = Y1
    1696:          t41 = Z1

    1698:          t12 = 0.0d0
    1699:          t22 = X2
    1700:          t32 = Y2
    1701:          t42 = Z2

    1703:          t13 =  c0*beta
    1704:          t23 = xnorm*beta + u0*(ubar0 + c0)
    1705:          t33 = ynorm*beta + v0*(ubar0 + c0)
    1706:          t43 = znorm*beta + w0*(ubar0 + c0)

    1708:          t14 = -c0*beta
    1709:          t24 = xnorm*beta + u0*(ubar0 - c0)
    1710:          t34 = ynorm*beta + v0*(ubar0 - c0)
    1711:          t44 = znorm*beta + w0*(ubar0 - c0)

    1713:          ti11 = -(u0*phi4 + v0*phi5 + w0*phi6)/beta
    1714:          ti21 = -(u0*phi7 + v0*phi8 + w0*phi9)/beta
    1715:          ti31 =  .5d0*(c0 - ubar0)/beta
    1716:          ti41 = -.5d0*(c0 + ubar0)/beta
    1717: !
    1718: ! Now, get the variables on the "inside"
    1719: !
    1720:          pi      = qnod(1,inode)
    1721:          ui      = qnod(2,inode)
    1722:          vi      = qnod(3,inode)
    1723:          wi      = qnod(4,inode)
    1724:          unorm   = xnorm*ui + ynorm*vi + znorm*wi
    1725: !
    1726: ! If ubar is negative, take the reference condition from outside
    1727: !
    1728: !
    1729:          if (unorm.gt.0.0d0)then
    1730:           pr = pi
    1731:           ur = ui
    1732:           vr = vi
    1733:           wr = wi
    1734:          else
    1735:           pr = p0
    1736:           ur = u0
    1737:           vr = v0
    1738:           wr = w0
    1739:          end if
    1740: !
    1741: ! Set rhs
    1742: !
    1743:          rhs1 = (ti11*pr + phi4*ur + phi5*vr + phi6*wr)/c20
    1744:          rhs2 = (ti21*pr + phi7*ur + phi8*vr + phi9*wr)/c20
    1745:          rhs3 = .5d0*(2.d0*ti31*pi + xnorm*ui + ynorm*vi + znorm*wi)/c20
    1746:          rhs4 = .5d0*(2.d0*ti41*p0 + xnorm*u0 + ynorm*v0 + znorm*w0)/c20
    1747: !
    1748: ! Now do matrix multiplication to get values on boundary
    1749: !
    1750:          pb =                       t13*rhs3 + t14*rhs4
    1751:          ub = t21*rhs1 + t22*rhs2 + t23*rhs3 + t24*rhs4
    1752:          vb = t31*rhs1 + t32*rhs2 + t33*rhs3 + t34*rhs4
    1753:          wb = t41*rhs1 + t42*rhs2 + t43*rhs3 + t44*rhs4

    1755:          ubar = xnorm*ub + ynorm*vb + znorm*wb
    1756:          flops = flops + 180.0
    1757: !
    1758:         if (inode.le.nnodes) then
    1759:           res(1,inode) = res(1,inode)+area*beta*ubar
    1760:           res(2,inode) = res(2,inode)+area*(ub*ubar + xnorm*pb)
    1761:           res(3,inode) = res(3,inode)+area*(vb*ubar + ynorm*pb)
    1762:           res(4,inode) = res(4,inode)+area*(wb*ubar + znorm*pb)
    1763:           flops = flops + 18.0
    1764:         endif
    1765: !
    1766:        enddo
    1767:        call PetscLogFlops(flops,ierr)
    1768: !
    1769: ! End of subroutine FLUX
    1770: !
    1771:       return
    1772:       end

    1774: !---------------------------------------------------------------
    1775: ! The following subroutines are from node3t.f in the original
    1776: ! code - D. K. Kaushik (1/17/97)
    1777: !---------------------------------------------------------------
    1778: !
    1779: !=============================== SUMGS ===============================72
    1780: !
    1781: ! Gets the weights for calculating gradients using least squares
    1782: !
    1783: !=====================================================================72
    1784:       subroutine SUMGS(nnodes,nedge,evec,xyz,rxy,irank,nvertices)

    1786:       implicit none
    1787: #include <petsc/finclude/petscsys.h>
    1788:       integer nnodes,nedge,irank,nvertices
    1789: #if defined(INTERLACING)
    1790:       PetscScalar xyz(3,nvertices)
    1791:       PetscScalar rxy(7,nnodes)
    1792:       integer evec(2,nedge)
    1793: #define x(i) xyz(1,i)
    1794: #define y(i) xyz(2,i)
    1795: #define z(i) xyz(3,i)
    1796: #define r11(i) rxy(1,i)
    1797: #define r12(i) rxy(2,i)
    1798: #define r13(i) rxy(3,i)
    1799: #define r22(i) rxy(4,i)
    1800: #define r23(i) rxy(5,i)
    1801: #define r33(i) rxy(6,i)
    1802: #define r44(i) rxy(7,i)
    1803: #undef eptr
    1804: #define eptr(j,i) evec(i,j)
    1805: #else
    1806:       PetscScalar xyz(nvertices,3)
    1807:       PetscScalar rxy(nnodes,7)
    1808:       integer evec(nedge,2)
    1809: #define x(i) xyz(i,1)
    1810: #define y(i) xyz(i,2)
    1811: #define z(i) xyz(i,3)
    1812: #define r11(i) rxy(i,1)
    1813: #define r12(i) rxy(i,2)
    1814: #define r13(i) rxy(i,3)
    1815: #define r22(i) rxy(i,4)
    1816: #define r23(i) rxy(i,5)
    1817: #define r33(i) rxy(i,6)
    1818: #define r44(i) rxy(i,7)
    1819: #define eptr(i,j) evec(i,j)
    1820: #endif
    1821:       integer i,n,node1,node2,ierr
    1822:       PetscScalar x1,y1,z1,x2,y2,z2,dx,dy,dz
    1823:       PetscScalar dist,weight,w2,w11,w22,w33
    1824:       PetscScalar r12r11,r13r11,r23r22,rmess

    1826: !
    1827: ! Initialize all the rij to 0.0
    1828: !
    1829:       do 1000 i = 1,nnodes
    1830:          r11(i) = 0.0
    1831:          r12(i) = 0.0
    1832:          r13(i) = 0.0
    1833:          r22(i) = 0.0
    1834:          r23(i) = 0.0
    1835:          r33(i) = 0.0
    1836:          r44(i) = 0.0
    1837:  1000 continue
    1838: !
    1839: ! Now loop over the edges and accumulate the r
    1840: !
    1841:       do 1020 n = 1, nedge
    1842: !
    1843:           node1 = eptr(n,1)
    1844:           node2 = eptr(n,2)
    1845: !
    1846:           x1 = x(node1)
    1847:           y1 = y(node1)
    1848:           z1 = z(node1)
    1849:           x2 = x(node2)
    1850:           y2 = y(node2)
    1851:           z2 = z(node2)
    1852: !
    1853:           dx = x2 - x1
    1854:           dy = y2 - y1
    1855:           dz = z2 - z1
    1856:           dist = sqrt(dx*dx + dy*dy + dz*dz)
    1857: !         weight = 1.0d0/dist
    1858:           weight = 1.0d0
    1859:           w2 = weight*weight
    1860: !
    1861:           if (node1 .le. nnodes) then
    1862:            r11(node1) = r11(node1) + (x2 - x1)*(x2 - x1)*w2
    1863:            r12(node1) = r12(node1) + (x2 - x1)*(y2 - y1)*w2
    1864:            r13(node1) = r13(node1) + (x2 - x1)*(z2 - z1)*w2
    1865:           endif
    1866: !
    1867:           if (node2 .le. nnodes) then
    1868:            r11(node2) = r11(node2) + (x1 - x2)*(x1 - x2)*w2
    1869:            r12(node2) = r12(node2) + (x1 - x2)*(y1 - y2)*w2
    1870:            r13(node2) = r13(node2) + (x1 - x2)*(z1 - z2)*w2
    1871:           endif
    1872:  1020 continue
    1873: !
    1874: !/*
    1875: ! Now calculate ||x(:)|| = r11 by taking the square root
    1876: ! Also divide r12 and r13 by ||x(:)||
    1877: !*/
    1878:       do 1030 i = 1,nnodes
    1879:         r11(i) = sqrt(r11(i))
    1880:         r12(i) = r12(i)/r11(i)
    1881:         r13(i) = r13(i)/r11(i)
    1882:  1030 continue
    1883: !
    1884: !/*
    1885: ! Now calculate r22 and r23
    1886: !*/
    1887:       do 1050 n = 1, nedge
    1888: !
    1889:           node1 = eptr(n,1)
    1890:           node2 = eptr(n,2)
    1891: !
    1892:           x1 = x(node1)
    1893:           y1 = y(node1)
    1894:           z1 = z(node1)
    1895:           x2 = x(node2)
    1896:           y2 = y(node2)
    1897:           z2 = z(node2)
    1898: !
    1899:           dx = x2 - x1
    1900:           dy = y2 - y1
    1901:           dz = z2 - z1
    1902:           dist = sqrt(dx*dx + dy*dy + dz*dz)
    1903: !         weight = 1.0/dist
    1904:           weight = 1.0d0
    1905:           dx = weight*dx
    1906:           dy = weight*dy
    1907:           dz = weight*dz
    1908:           w2 = weight*weight
    1909: !
    1910:           if (node1 .le. nnodes) then
    1911:            r22(node1) = r22(node1) +                                         &
    1912:      &                  (dy-dx*r12(node1)/r11(node1))**2
    1913:            r23(node1) = r23(node1) + dz*                                     &
    1914:      &                  (dy-dx*r12(node1)/r11(node1))
    1915:           endif
    1916: !
    1917:           if (node2 .le. nnodes) then
    1918:            r22(node2) = r22(node2) +                                         &
    1919:      &                  (-dy + dx*r12(node2)/r11(node2))**2
    1920:            r23(node2) = r23(node2) - dz*                                     &
    1921:      &                  (-dy+dx*r12(node2)/r11(node2))
    1922:           endif
    1923:  1050 continue
    1924: !
    1925: !/*
    1926: ! Now finish getting r22 and r23
    1927: !*/
    1928:       do 1060 i = 1,nnodes
    1929:         r22(i) = sqrt(r22(i))
    1930:         r23(i) = r23(i)/r22(i)
    1931:  1060 continue
    1932: !
    1933: !/*
    1934: ! Now all we have to do is get r33
    1935: !*/
    1936:       do 1080 n = 1, nedge
    1937: !
    1938:           node1 = eptr(n,1)
    1939:           node2 = eptr(n,2)
    1940: !
    1941:           x1 = x(node1)
    1942:           y1 = y(node1)
    1943:           z1 = z(node1)
    1944:           x2 = x(node2)
    1945:           y2 = y(node2)
    1946:           z2 = z(node2)
    1947: !
    1948:           dx = x2 - x1
    1949:           dy = y2 - y1
    1950:           dz = z2 - z1
    1951:           dist = sqrt(dx*dx + dy*dy + dz*dz)
    1952: !         weight = 1.0/dist
    1953:           weight = 1.0d0
    1954:           dx = weight*dx
    1955:           dy = weight*dy
    1956:           dz = weight*dz
    1957:           w2 = weight*weight
    1958: !
    1959:           if (node1 .le. nnodes) then
    1960:            r33(node1) = r33(node1) +                                        &
    1961:      &                  (dz-dx*r13(node1)/r11(node1)-                       &
    1962:      &                  r23(node1)/r22(node1)*                              &
    1963:      &                  (dy - dx*r12(node1)/                                &
    1964:      &                   r11(node1)))**2
    1965:           endif
    1966: !
    1967:           if (node2 .le. nnodes) then
    1968:            r33(node2) = r33(node2) +                                        &
    1969:      &                  (-dz+dx*r13(node2)/r11(node2)-                      &
    1970:      &                  r23(node2)/r22(node2)*                              &
    1971:      &                  (-dy + dx*r12(node2)/                               &
    1972:      &                  r11(node2)))**2
    1973:           endif
    1974: !
    1975:  1080 continue
    1976: !
    1977: !/*
    1978: ! Now just get the magnitude of r33
    1979: !*/
    1980:       do 1090 i = 1,nnodes
    1981:         r33(i) = sqrt(r33(i))
    1982:  1090 continue
    1983: !
    1984: !/*
    1985: ! Added by Dinesh Kaushik (6/27/97)
    1986: ! The following addition changes the meaning of r11 .. r33. r44
    1987: ! is the new parameter introduced by me. The new definitions
    1988: ! are taken from LSTGS (where these parameters
    1989: ! are used finally).
    1990: ! r11->w11
    1991: ! r22->w22
    1992: ! r33->w33
    1993: ! r12->r12r11
    1994: ! r13->r13r11
    1995: ! r23->r23r22
    1996: ! r44->rmess */

    1998:       do i = 1,nnodes
    1999:            w11 = 1.d0/(r11(i)*r11(i))
    2000:            w22 = 1.d0/(r22(i)*r22(i))
    2001:            w33 = 1.d0/(r33(i)*r33(i))
    2002:            r12r11 = r12(i)/r11(i)
    2003:            r13r11 = r13(i)/r11(i)
    2004:            r23r22 = r23(i)/r22(i)
    2005:            rmess=(r12(i)*r23(i)-                                             &
    2006:      &            r13(i)*r22(i))/                                            &
    2007:      &           (r11(i)*r22(i)*                                             &
    2008:      &            r33(i)*r33(i))

    2010:            r11(i) = w11
    2011:            r22(i) = w22
    2012:            r33(i) = w33
    2013:            r12(i) = r12r11
    2014:            r13(i) = r13r11
    2015:            r23(i) = r23r22
    2016:            r44(i) = rmess
    2017:       enddo
    2018: !
    2019: ! Finished with SUMGS
    2020: !
    2021:       return
    2022:       end

    2025: !================================= LSTGS =============================72
    2026: !
    2027: ! Calculates the Gradients at the nodes using weighted least squares
    2028: ! This subroutine solves using Gram-Schmidt
    2029: !
    2030: !=====================================================================72
    2031:       subroutine LSTGS(nnodes,nedge,evec,                                &
    2032:      &                 qvec,grad,xyz,                                    &
    2033:      &                 rxy,irank,nvertices)
    2034: !
    2035:       implicit none
    2036: #include <petsc/finclude/petscsys.h>
    2037:       integer nnodes,nedge,irank,nvertices
    2038: !
    2039: #if defined(INTERLACING)
    2040:       PetscScalar qvec(4,nvertices)
    2041:       PetscScalar xyz(3,nvertices)
    2042:       PetscScalar rxy(7,nnodes)
    2043:       PetscScalar grad(3,4,nnodes)
    2044:       integer evec(2,nedge)
    2045: #define qnod(i,j) qvec(i,j)
    2046: #define gradx(i,j) grad(1,i,j)
    2047: #define grady(i,j) grad(2,i,j)
    2048: #define gradz(i,j) grad(3,i,j)
    2049: #define x(i) xyz(1,i)
    2050: #define y(i) xyz(2,i)
    2051: #define z(i) xyz(3,i)
    2052: #define r11(i) rxy(1,i)
    2053: #define r12(i) rxy(2,i)
    2054: #define r13(i) rxy(3,i)
    2055: #define r22(i) rxy(4,i)
    2056: #define r23(i) rxy(5,i)
    2057: #define r33(i) rxy(6,i)
    2058: #define r44(i) rxy(7,i)
    2059: #define eptr(j,i) evec(i,j)
    2060: #else
    2061:       PetscScalar qvec(nvertices,4)
    2062:       PetscScalar xyz(nvertices,3)
    2063:       PetscScalar rxy(nnodes,7)
    2064:       PetscScalar grad(nnodes,4,3)
    2065:       integer evec(nedge,2)
    2066: #define qnod(i,j) qvec(j,i)
    2067: #define gradx(i,j) grad(j,i,1)
    2068: #define grady(i,j) grad(j,i,2)
    2069: #define gradz(i,j) grad(j,i,3)
    2070: #define x(i) xyz(i,1)
    2071: #define y(i) xyz(i,2)
    2072: #define z(i) xyz(i,3)
    2073: #define r11(i) rxy(i,1)
    2074: #define r12(i) rxy(i,2)
    2075: #define r13(i) rxy(i,3)
    2076: #define r22(i) rxy(i,4)
    2077: #define r23(i) rxy(i,5)
    2078: #define r33(i) rxy(i,6)
    2079: #define r44(i) rxy(i,7)
    2080: #define eptr(i,j) evec(i,j)
    2081: #endif
    2082:       PetscScalar  title(20),beta,alpha
    2083:       PetscScalar  Re,dt,tot,res0,resc
    2084:       integer i,ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    2085:       PetscScalar  cfl1,cfl2
    2086:       integer nsmoth,iflim,itran,nbtran,jupdate,                          &
    2087:      &        nstage,ncyct,iramp,nitfo
    2088:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                   &
    2089:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    2090:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,           &
    2091:      &             nstage,ncyct,iramp,nitfo
    2092:       integer n,node1,node2,ierr
    2093:       PetscLogDouble flops
    2094:       PetscScalar  dx1,dy1,dz1,dx2,dy2,dz2,                                    &
    2095:      &        dq1,dq2,dq3,dq4,                                            &
    2096:      &        weight,w11,r12r11,r13r11,w22,r23r22,w33,rmess,              &
    2097:      &        coef1,coef2,termx,termy,termz
    2098: !     logging variables
    2099: !     integer flag
    2100: !     integer grad_event,node1_event,node2_event
    2101: !     character * 16 grad_label, node1_label, node2_label
    2102: !     data flag/-1/,grad_label/'GRAD            '/
    2103: !     data node1_label/'NODE1           '/
    2104: !     data node2_label/'NODE2           '/
    2105: !     save grad_event, grad_label, flag
    2106: !     save node1_event,node2_event,node1_label, node2_label

    2108: !     if (flag .eq. -1) then
    2109: !        call PetscLogEventRegister(grad_event,grad_label,ierr)
    2110: !        call PetscLogEventRegister(node1_event,node1_label,ierr)
    2111: !        call PetscLogEventRegister(node2_event,node2_label,ierr)
    2112: !        flag = 1
    2113: !     endif
    2114: !     call PetscLogEventBegin(grad_event,0,0,0,0,ierr)

    2116:       flops = 0.0
    2117: ! For checking out the code input a linear distribution
    2118: !
    2119: !     write(6,700)nnodes,nedge
    2120: ! 700 format(1h ,'nnodes=',i5,' nedge=',i5)
    2121: !     do 1001 i = 1,nnodes
    2122: !     write(6,800)i,x(i),y(i),z(i)
    2123: ! 800 format(1h ,'i x y z=',i5,3(f10.5,1x))
    2124: !       qnod(1,i) = 1.0*x(i) +  2.0*y(i) + 3.0*z(i)
    2125: !       qnod(2,i) = 3.0*x(i) +  4.0*y(i) + 6.0*z(i)
    2126: !       qnod(3,i) = 5.0*x(i) +  6.0*y(i) + 9.0*z(i)
    2127: !       qnod(4,i) = 7.0*x(i) +  8.0*y(i) + 12.0*z(i)
    2128: !       qnod(i,5) = 9.0*x(i) + 10.0*y(i) + 15.0*z(i)
    2129: !1001 continue
    2130: !
    2131: ! Zero out the gradients
    2132: !
    2133:       do 1000 i = 1,nnodes
    2134:          gradx(1,i) = 0.0
    2135:          grady(1,i) = 0.0
    2136:          gradz(1,i) = 0.0

    2138:          gradx(2,i) = 0.0
    2139:          grady(2,i) = 0.0
    2140:          gradz(2,i) = 0.0

    2142:          gradx(3,i) = 0.0
    2143:          grady(3,i) = 0.0
    2144:          gradz(3,i) = 0.0

    2146:          gradx(4,i) = 0.0
    2147:          grady(4,i) = 0.0
    2148:          gradz(4,i) = 0.0
    2149:  1000 continue

    2151: ! If second order, loop over all the faces accumulate sums
    2152: !
    2153: !     nitfo = 0
    2154: !     if (ntt.gt.nitfo.or.ivisc.gt.0)then
    2155: !     if (1 .gt. 0) then
    2156:        do 1020 n = 1, nedge
    2157:          node1 = eptr(n,1)
    2158:          node2 = eptr(n,2)
    2159: !        if ((node1 .le. nnodes).or.(node2 .le. nnodes)) then
    2160:            dx1 = x(node2) - x(node1)
    2161:            dy1 = y(node2) - y(node1)
    2162:            dz1 = z(node2) - z(node1)
    2163: !
    2164: !          dist = sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1)
    2165: !          weight = 1.0/dist
    2166:            weight = 1.0d0
    2167: !          w2 = weight*weight
    2168: !
    2169:            dx1 = weight*dx1
    2170:            dy1 = weight*dy1
    2171:            dz1 = weight*dz1

    2173:            flops = flops + 6.0
    2174: !
    2175: !         call PetscLogEventBegin(node1_event,0,0,0,0,ierr)
    2176:           if (node1 .le. nnodes) then
    2177:            dq1 = weight*(qnod(1,node2) - qnod(1,node1))
    2178:            dq2 = weight*(qnod(2,node2) - qnod(2,node1))
    2179:            dq3 = weight*(qnod(3,node2) - qnod(3,node1))
    2180:            dq4 = weight*(qnod(4,node2) - qnod(4,node1))
    2181: !
    2182: !          w11 = 1./(r11(node1)*r11(node1))
    2183: !          w22 = 1./(r22(node1)*r22(node1))
    2184: !          w33 = 1./(r33(node1)*r33(node1))
    2185: !          r12r11 = r12(node1)/r11(node1)
    2186: !          r13r11 = r13(node1)/r11(node1)
    2187: !          r23r22 = r23(node1)/r22(node1)
    2188: !          rmess  = (r12(node1)*r23(node1) - r13(node1)*r22(node1))/
    2189: !    1              (r11(node1)*r22(node1)*r33(node1)*r33(node1))
    2190: !
    2191:            w11 = r11(node1)
    2192:            r12r11 = r12(node1)
    2193:            r13r11 = r13(node1)
    2194:            w22 = r22(node1)
    2195:            r23r22 = r23(node1)
    2196:            w33 = r33(node1)
    2197:            rmess  = r44(node1)
    2198: !
    2199:            coef1  = dy1 - dx1*r12r11
    2200:            coef2  = dz1 - dx1*r13r11 - r23r22*coef1
    2201:            termx = dx1*w11 - w22*r12r11*coef1 + rmess*coef2
    2202:            termy = w22*coef1 - r23r22*w33*coef2
    2203:            termz = w33*coef2
    2204: !
    2205:            gradx(1,node1) = gradx(1,node1) + termx*dq1
    2206:            grady(1,node1) = grady(1,node1) + termy*dq1
    2207:            gradz(1,node1) = gradz(1,node1) + termz*dq1
    2208: !
    2209:            gradx(2,node1) = gradx(2,node1) + termx*dq2
    2210:            grady(2,node1) = grady(2,node1) + termy*dq2
    2211:            gradz(2,node1) = gradz(2,node1) + termz*dq2
    2212: !
    2213:            gradx(3,node1) = gradx(3,node1) + termx*dq3
    2214:            grady(3,node1) = grady(3,node1) + termy*dq3
    2215:            gradz(3,node1) = gradz(3,node1) + termz*dq3
    2216: !
    2217:            gradx(4,node1) = gradx(4,node1) + termx*dq4
    2218:            grady(4,node1) = grady(4,node1) + termy*dq4
    2219:            gradz(4,node1) = gradz(4,node1) + termz*dq4
    2220: !
    2221:            flops = flops + 49.0
    2222:           endif
    2223: !         call PetscLogEventEnd(node1_event,0,0,0,0,ierr)
    2224: !
    2225: ! Now do the other node
    2226: !
    2227: !         call PetscLogEventBegin(node2_event,0,0,0,0,ierr)
    2228:           if (node2 .le. nnodes) then
    2229:            dx2 = -dx1
    2230:            dy2 = -dy1
    2231:            dz2 = -dz1
    2232: !
    2233:            dq1 = weight*(qnod(1,node1) - qnod(1,node2))
    2234:            dq2 = weight*(qnod(2,node1) - qnod(2,node2))
    2235:            dq3 = weight*(qnod(3,node1) - qnod(3,node2))
    2236:            dq4 = weight*(qnod(4,node1) - qnod(4,node2))
    2237: !
    2238: !          w11 = 1./(r11(node2)*r11(node2))
    2239: !          w22 = 1./(r22(node2)*r22(node2))
    2240: !          w33 = 1./(r33(node2)*r33(node2))
    2241: !          r12r11 = r12(node2)/r11(node2)
    2242: !          r13r11 = r13(node2)/r11(node2)
    2243: !          r23r22 = r23(node2)/r22(node2)
    2244: !          rmess  = (r12(node2)*r23(node2) - r13(node2)*r22(node2))/
    2245: !    &              (r11(node2)*r22(node2)*r33(node2)*r33(node2))

    2247:            w11 = r11(node2)
    2248:            r12r11 = r12(node2)
    2249:            r13r11 = r13(node2)
    2250:            w22 = r22(node2)
    2251:            r23r22 = r23(node2)
    2252:            w33 = r33(node2)
    2253:            rmess  = r44(node2)
    2254: !
    2255:            coef1  = dy2 - dx2*r12r11
    2256:            coef2  = dz2 - dx2*r13r11 - r23r22*coef1
    2257:            termx = dx2*w11 - w22*r12r11*coef1 + rmess*coef2
    2258:            termy = w22*coef1 - r23r22*w33*coef2
    2259:            termz = w33*coef2
    2260: !
    2261:            gradx(1,node2) = gradx(1,node2) + termx*dq1
    2262:            grady(1,node2) = grady(1,node2) + termy*dq1
    2263:            gradz(1,node2) = gradz(1,node2) + termz*dq1
    2264: !
    2265:            gradx(2,node2) = gradx(2,node2) + termx*dq2
    2266:            grady(2,node2) = grady(2,node2) + termy*dq2
    2267:            gradz(2,node2) = gradz(2,node2) + termz*dq2
    2268: !
    2269:            gradx(3,node2) = gradx(3,node2) + termx*dq3
    2270:            grady(3,node2) = grady(3,node2) + termy*dq3
    2271:            gradz(3,node2) = gradz(3,node2) + termz*dq3
    2272: !
    2273:            gradx(4,node2) = gradx(4,node2) + termx*dq4
    2274:            grady(4,node2) = grady(4,node2) + termy*dq4
    2275:            gradz(4,node2) = gradz(4,node2) + termz*dq4
    2276: !
    2277:            flops = flops + 52.0
    2278:           endif
    2279: !         call PetscLogEventEnd(node2_event,0,0,0,0,ierr)
    2280: !       endif
    2281: !
    2282:  1020  continue
    2283: !     end if
    2284:       call PetscLogFlops(flops,ierr)
    2285: !     call PetscLogEventEnd(grad_event,0,0,0,0,ierr)
    2286: !
    2287: ! End of LSTGS
    2288: !
    2289:       return
    2290:       end

    2293: !=================================== GETRES ==========================72
    2294: !
    2295: ! Calculates the residual
    2296: ! Last Modified - D. K. Kaushik 1/23/97
    2297: ! I have eliminated the input variables which were not needed -
    2298: ! dq, A, B, iupdate
    2299: !
    2300: !=====================================================================72
    2301: #if defined(_OPENMP)
    2302: #if defined(HAVE_EDGE_COLORING)
    2303:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
    2304:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
    2305:      &                  ileast,isnode,ivnode,ifnode,                        &
    2306:      &                  nnfacet,f2ntn,nnbound,                              &
    2307:      &                  nvfacet,f2ntv,nvbound,                              &
    2308:      &                  nffacet,f2ntf,nfbound,                              &
    2309:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
    2310:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
    2311:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
    2312:      &                  max_threads,                                        &
    2313:      &                  ncolor,ncount,                                      &
    2314:      &                  LocalTS,irank, nvertices)
    2315: #elif defined(HAVE_REDUNDANT_WORK)
    2316:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
    2317:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
    2318:      &                  ileast,isnode,ivnode,ifnode,                        &
    2319:      &                  nnfacet,f2ntn,nnbound,                              &
    2320:      &                  nvfacet,f2ntv,nvbound,                              &
    2321:      &                  nffacet,f2ntf,nfbound,                              &
    2322:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
    2323:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
    2324:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
    2325:      &                  max_threads,                                        &
    2326:      &                  resd,                                               &
    2327:      &                  LocalTS,irank, nvertices)
    2328: #else
    2329:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
    2330:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
    2331:      &                  ileast,isnode,ivnode,ifnode,                        &
    2332:      &                  nnfacet,f2ntn,nnbound,                              &
    2333:      &                  nvfacet,f2ntv,nvbound,                              &
    2334:      &                  nffacet,f2ntf,nfbound,                              &
    2335:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
    2336:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
    2337:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
    2338:      &                  max_threads,                                        &
    2339:      &                  nedgeAllThr,part_thr,nedge_thr,                     &
    2340:      &                  edge_thr,xyzn_thr,                                  &
    2341:      &                  LocalTS,irank, nvertices)
    2342: #endif
    2343: #else
    2344:       subroutine GETRES(nnodes,ncell,nedge,nsface,nvface,nfface,nbface,     &
    2345:      &                  nsnode,nvnode,nfnode,isface,ivface,ifface,          &
    2346:      &                  ileast,isnode,ivnode,ifnode,                        &
    2347:      &                  nnfacet,f2ntn,nnbound,                              &
    2348:      &                  nvfacet,f2ntv,nvbound,                              &
    2349:      &                  nffacet,f2ntf,nfbound,                              &
    2350:      &                  evec,sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,           &
    2351:      &                  xyzn,qvec,cdt,xyz,area,grad,resvec,turbre,          &
    2352:      &                  slen,c2n,c2e,us,vs,as,phi,amut,ires,                &
    2353:      &                  LocalTS,irank, nvertices)
    2354: #endif
    2355: !
    2356:       implicit none
    2357: #include <petsc/finclude/petscsys.h>
    2358: !
    2359:       integer nnodes, ncell, nedge,                                         &
    2360:      &        nbface,ileast,ires,                                           &
    2361:      &        nsface, nvface, nfface,                                       &
    2362:      &        nsnode, nvnode, nfnode,                                       &
    2363:      &        nnfacet,nvfacet,nffacet,                                      &
    2364:      &        nnbound,nvbound,nfbound,                                      &
    2365:      &        LocalTS,irank,nvertices
    2366:       integer evec(nedge,2)
    2367:       integer isface(1),ivface(1),ifface(1)
    2368:       integer isnode(1),ivnode(1),ifnode(1)
    2369:       integer c2n(ncell,4),c2e(ncell,6)
    2370:       integer f2ntn(nnfacet,4)
    2371:       integer f2ntv(nvfacet,4)
    2372:       integer f2ntf(nffacet,4)
    2373: !
    2374:       PetscScalar us(nbface,3,4),vs(nbface,3,4)
    2375:       PetscScalar as(nbface,3,4)
    2376:       PetscScalar sxn(1),syn(1),szn(1)
    2377:       PetscScalar vxn(1),vyn(1),vzn(1)
    2378:       PetscScalar fxn(1),fyn(1),fzn(1)
    2379:       PetscScalar xyz(nvertices,3),area(nvertices)
    2380:       PetscScalar xyzn(nedge,4)
    2381:       PetscScalar turbre(1),slen(1)
    2382:       PetscScalar qvec(nvertices,4)
    2383:       PetscScalar phi(nvertices,4)
    2384:       PetscScalar cdt(nvertices)
    2385:       PetscScalar grad(3,4,nvertices)
    2386: !     real dq(nnodes,4)
    2387: !     real r11(nvertices),r12(nvertices),r13(nvertices)
    2388: !     real r22(nvertices),r23(nvertices),r33(nvertices)
    2389:       PetscScalar amut(nnodes)
    2390: !
    2391:       integer i
    2392:       PetscScalar  title(20),beta,alpha
    2393:       PetscScalar  Re,dt,tot,res0,resc
    2394:       integer ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    2395:       PetscScalar  cfl1,cfl2
    2396:       integer nsmoth,iflim,itran,nbtran,jupdate,                             &
    2397:      &        nstage,ncyct,iramp,nitfo
    2398:       PetscScalar  gtol
    2399:       integer icycle,nsrch,ilu0,ifcn
    2400:       PetscScalar  rms(1001),clw(1001),cdw(1001),                                 &
    2401:      &        cmw(1001),xres(1001)
    2402:       common/info/title,beta,alpha,Re,dt,tot,res0,resc,                      &
    2403:      &            ntt,mseq,ivisc,irest,icyc,ihane,ntturb
    2404:       common/runge/cfl1,cfl2,nsmoth,iflim,itran,nbtran,jupdate,              &
    2405:      &             nstage,ncyct,iramp,nitfo
    2406:       common/gmcom/gtol,icycle,nsrch,ilu0,ifcn
    2407:       common/history/rms,clw,cdw,cmw,xres
    2408: !     logging variables
    2409: !     integer flux_event, flag, ierr
    2410: !     integer delta2_event
    2411: !     character * 16 flux_label
    2412: !     character * 16 delta2_label
    2413: !     data flag/-1/,flux_label/'FluxEvaluation  '/
    2414: !     data delta2_label/'DELTA2          '/
    2415: !     save flux_event, flag, flux_label
    2416: !     save delta2_event,delta2_label
    2417: #if defined(INTERLACING)
    2418:        PetscScalar resvec(4,nnodes)
    2419: #undef res
    2420: #define res(i,j) resvec(i,j)
    2421: #else
    2422:        PetscScalar resvec(nnodes,4)
    2423: #define res(i,j) resvec(j,i)
    2424: #endif
    2425: #if defined(_OPENMP)
    2426:        integer max_threads
    2427: #if defined(HAVE_EDGE_COLORING)
    2428:        integer ncolor,ncount(1)
    2429: #elif defined(HAVE_REDUNDANT_WORK)
    2430:        PetscScalar resd(4,nnodes)
    2431: #else
    2432:        integer nedgeAllThr,part_thr(nnodes),nedge_thr(0:max_threads),         &
    2433:      &         edge_thr(2,nedgeAllThr)
    2434:        PetscScalar xyzn_thr(4,nedgeAllThr)
    2435: #endif
    2436: #endif
    2437: !

    2439: !     if (flag .eq. -1) then
    2440: !        call PetscLogEventRegister(delta2_event,delta2_label,ierr)
    2441: !        call PetscLogEventRegister(flux_event,flux_label,ierr)
    2442: !        flag = 1
    2443: !     endif
    2444: !
    2445: ! Calculate the time step
    2446: !
    2447:       if (ires.eq.1) goto 888
    2448: !     call PetscLogEventBegin(delta2_event,0,0,0,0,ierr)
    2449:       call DELTAT2(nnodes,nedge,qvec,cdt,                                   &
    2450:      &             xyz,area,xyzn,evec,                                      &
    2451:      &             sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,                     &
    2452:      &             nsnode,nvnode,nfnode,isnode,ivnode,ifnode,               &
    2453:      &             LocalTS,irank,nvertices)
    2454: !     call PetscLogEventEnd(delta2_event,0,0,0,0,ierr)
    2455:   888 continue
    2456: !
    2457: !/*
    2458: !   Calculate the gradients
    2459: !   ----Kyle seems to recommend only LSTGS for gradients,
    2460: !   so I have commented the GETGRAD call - DKK (1/17/97)
    2461: !
    2462: !     if (ileast.eq.0) then
    2463: !        call GETGRAD(nnodes,ncell,nedge,nsface,nvface,nfface,
    2464: !    &                isface,ivface,ifface,eptr,ncolor,ncount,
    2465: !    &                qnod,gradx,grady,x,y,
    2466: !    &                area,wx,wy,xn,yn,rl)
    2467: !
    2468: !     else if (ileast.eq.4) then
    2469: !     if (ileast.eq.4) then
    2470: !        call LSTGS(nnodes,nedge,eptr,
    2471: !    &              qnod,gradx,grady,gradz,xyz,
    2472: !    &              r11,r12,r13,r22,r23,r33,irank,nvertices)
    2473: !     end if
    2474: !
    2475: ! zero out residuals (viscous residuals are zeroed in vfluxnew)
    2476: !*/
    2477:       do 1002 i = 1,nnodes
    2478:           res(1,i)=0.0d0
    2479:           res(2,i)=0.0d0
    2480:           res(3,i)=0.0d0
    2481:           res(4,i)=0.0d0

    2483:           phi(i,1)=1.0d0
    2484:           phi(i,2)=1.0d0
    2485:           phi(i,3)=1.0d0
    2486:           phi(i,4)=1.0d0
    2487:  1002   continue
    2488: !
    2489: !/*
    2490: !
    2491: ! If not doing Newton-Krylov and iflim=1 call the Flux Limiter
    2492: !
    2493: !     if (iflim.eq.1.and.ifcn.ne.1)then
    2494: !      call TIMLIM(nnodes,nedge,qnod,res,dq,phi,ncolor,ncount,
    2495: !    &             gradx,grady,gradz,x,y,z,eptr)
    2496: !
    2497: ! If we used the limiter we need to zero out the residual again
    2498: ! since we used it for scratch space
    2499: !
    2500: !     do 1003 i = 1,nnodes
    2501: !         res(1,i)=0.0
    2502: !         res(2,i)=0.0
    2503: !         res(3,i)=0.0
    2504: !         res(4,i)=0.0
    2505: !1003   continue
    2506: !
    2507: !     end if
    2508: !*/
    2509: !   Split the fluxes and perform the flux balance
    2510: !
    2511: !       call PetscLogEventBegin(flux_event,0,0,0,0,ierr)

    2513: #if defined(_OPENMP)
    2514: #if defined(HAVE_EDGE_COLORING)
    2515:       call FLUX(nnodes, ncell, nedge,                                     &
    2516:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
    2517:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    2518:      &          nnfacet,f2ntn,nnbound,                                    &
    2519:      &          nvfacet,f2ntv,nvbound,                                    &
    2520:      &          nffacet,f2ntf,nfbound,                                    &
    2521:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
    2522:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    2523:      &          max_threads,                                              &
    2524:      &          ncolor,ncount,                                            &
    2525:      &          irank,nvertices)
    2526: #elif defined(HAVE_REDUNDANT_WORK)
    2527:       call FLUX(nnodes, ncell, nedge,                                     &
    2528:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
    2529:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    2530:      &          nnfacet,f2ntn,nnbound,                                    &
    2531:      &          nvfacet,f2ntv,nvbound,                                    &
    2532:      &          nffacet,f2ntf,nfbound,                                    &
    2533:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
    2534:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    2535:      &          max_threads,                                              &
    2536:      &          resd,                                                     &
    2537:      &          irank,nvertices)
    2538: #else
    2539:       call FLUX(nnodes, ncell, nedge,                                     &
    2540:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
    2541:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    2542:      &          nnfacet,f2ntn,nnbound,                                    &
    2543:      &          nvfacet,f2ntv,nvbound,                                    &
    2544:      &          nffacet,f2ntf,nfbound,                                    &
    2545:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
    2546:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    2547:      &          max_threads,                                              &
    2548:      &          nedgeAllThr,part_thr,nedge_thr,edge_thr,xyzn_thr,         &
    2549:      &          irank,nvertices)
    2550: #endif
    2551: #else
    2552:       call FLUX(nnodes, ncell, nedge,                                     &
    2553:      &          nsface, nvface, nfface, isface, ivface, ifface,           &
    2554:      &          nsnode, nvnode, nfnode, isnode, ivnode, ifnode,           &
    2555:      &          nnfacet,f2ntn,nnbound,                                    &
    2556:      &          nvfacet,f2ntv,nvbound,                                    &
    2557:      &          nffacet,f2ntf,nfbound,                                    &
    2558:      &          grad,evec,qvec,xyz,resvec,xyzn,                           &
    2559:      &          sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,phi,                  &
    2560:      &          irank,nvertices)
    2561: #endif
    2562: !          call PetscLogEventEnd(flux_event,0,0,0,0,ierr)
    2563: !/*
    2564: ! calculate viscous fluxes
    2565: !
    2566: !     if (ivisc.gt.0) then
    2567: !           call VISRHS (nnodes,ncell,nedge,
    2568: !           call EDGEVIS(nnodes,ncell,nedge,
    2569: !    &                  nsnode,nvnode,nfnode,isnode,ivnode,ifnode,
    2570: !    &                  nsface,nvface,nfface,isface,ivface,ifface,
    2571: !    &                  nnfacet,f2ntn,nnbound,ncolorn,countn,
    2572: !    &                  nvfacet,f2ntv,nvbound,ncolorv,countv,
    2573: !    &                  nffacet,f2ntf,nfbound,ncolorf,countf,
    2574: !    &                  nccolor,nccount,
    2575: !    &                  eptr,c2n,c2e,
    2576: !    &                  sxn,syn,szn,vxn,vyn,vzn,fxn,fyn,fzn,
    2577: !    &                  x,y,z,gradx,grady,gradz,
    2578: !    &                  qnod,amut,res,phi)
    2579: !     end if
    2580: !*/
    2581: ! End of subroutine GETRES
    2582: !
    2583:       return
    2584:       end

    2586: !---------------------------------------------------------------
    2587: ! The following subroutine is from node4t.f in the original
    2588: ! code - D. K. Kaushik (1/17/97)
    2589: !---------------------------------------------------------------
    2590: !
    2591: !=============================================================================
    2592: !
    2593: !  Opens files for I/O
    2594: !
    2595: !=============================================================================
    2596:       SUBROUTINE OPENM(irank)
    2597: !
    2598: !  TAPE7  -- input:  mach number, angle of attack etc..
    2599: !  TAPE9  -- input:  reads restart file
    2600: !  TAPE10 -- output: residual history
    2601: !  TAPE11 -- output: writes restart file
    2602: !  TAPE12 -- output: writes residual and lift for plotting
    2603: !  TAPE13 -- output: writes flowfield for contour plotting
    2604: !
    2605:       implicit none
    2606:       integer  irank
    2607: !     OPEN(UNIT= 7,FILE='home/petsc/datafiles/fun3dgrid/ginput.faces',            &
    2608: !    &form='formatted',STATUS='OLD')

    2610: !     OPEN(UNIT=9,FILE='framer.bin',
    2611: !    &form='unformatted',STATUS='old')

    2613:       if (irank .eq. 0) OPEN(UNIT= 10,FILE='frame.out3',                          &
    2614:      &form='formatted',STATUS='unknown')

    2616: !     OPEN(UNIT= 11,FILE='frame.bin',
    2617: !    +form='unformatted',STATUS='unknown')

    2619: !     OPEN(UNIT= 12,FILE='frame.plt',
    2620: !    +form='formatted',STATUS='unknown')

    2622: !     OPEN(UNIT= 13,FILE='frame.tec',
    2623: !    +form='formatted',STATUS='unknown')

    2625: !     OPEN(UNIT= 14,FILE='frame.fast.g',
    2626: !    +form='unformatted',STATUS='unknown')

    2628: !     OPEN(UNIT= 15,FILE='frame.fast.q',
    2629: !    +form='unformatted',STATUS='unknown')

    2631:       return
    2632:       end
    2633: !
    2634: !
    2635: !===================================================================
    2636: !
    2637: ! Get the IA, JA, and IAU arrays
    2638: !
    2639: !===================================================================
    2640:       subroutine GETIA(nnodes,nedge,evec,ia,ideg,irank)
    2641:       implicit none
    2642: #include <petsc/finclude/petscsys.h>
    2643:       integer nnodes,nedge,irank
    2644:       integer ia(1),ideg(1)
    2645: #if defined(INTERLACING)
    2646:        integer evec(2,nedge)
    2647: #define eptr(j,i) evec(i,j)
    2648: #else
    2649:        integer evec(nedge,2)
    2650: #define eptr(i,j) evec(i,j)
    2651: #endif
    2652:        integer i,node1,node2
    2653: !
    2654: ! First get the degree of each node using ideg as a dummy array
    2655: !
    2656:       do 1000 i = 1,nnodes
    2657:         ideg(i) = 0
    2658:  1000 continue
    2659: !
    2660:       do 1010 i = 1,nedge
    2661:         node1 = eptr(i,1)
    2662:         node2 = eptr(i,2)
    2663: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
    2664:         ideg(node1) = ideg(node1) + 1
    2665:         ideg(node2) = ideg(node2) + 1
    2666: #else
    2667:         if (node1 .le. nnodes) ideg(node1) = ideg(node1) + 1
    2668:         if (node2 .le. nnodes) ideg(node2) = ideg(node2) + 1
    2669: #endif
    2670:  1010 continue
    2671: !
    2672: ! Now we can fill the ia array fairly easily
    2673: !
    2674:       ia(1) = 1
    2675:       do 1020 i = 1,nnodes
    2676:         ia(i+1) = ia(i) + ideg(i) + 1
    2677: !       write(9,100)i,ideg(i)
    2678: ! 100   format(1h ,'deg(',i6,')=',i6)
    2679:  1020 continue
    2680: !
    2681:       return
    2682:       end
    2683: !
    2684: !===================================================================
    2685: !
    2686: ! Get the IA, JA, and IAU arrays
    2687: !
    2688: !===================================================================
    2689:       subroutine GETJA(nnodes,nedge,evec,ia,ja,iwork,irank)
    2690:       implicit none
    2691: #include <petsc/finclude/petscsys.h>
    2692:       integer nnodes,nedge,irank
    2693:       integer ia(1),ja(1),iwork(1)
    2694: #if defined(INTERLACING)
    2695:        integer evec(2,nedge)
    2696: #define eptr(j,i) evec(i,j)
    2697: #else
    2698:        integer evec(nedge,2)
    2699: #define eptr(i,j) evec(i,j)
    2700: #endif
    2701:        integer i,index,node1,node2,index1,index2,                          &
    2702:      &         istart,iend
    2703: !     open(unit=90,file='map.out',status='UNKNOWN')
    2704: !
    2705: ! Now we need to get the JA array
    2706: ! First fill the diagonal places
    2707: !
    2708:       do 1040 i = 1,nnodes
    2709:         index = ia(i)
    2710:         ja(index) = i
    2711:         iwork(i) = 1
    2712:  1040 continue
    2713: !
    2714:       do 1030 i = 1,nedge
    2715:         node1 = eptr(i,1)
    2716:         node2 = eptr(i,2)
    2717: !
    2718: #if defined(_OPENMP) && !defined(HAVE_REDUNDANT_WORK) && !defined(HAVE_EDGE_COLORING)
    2719:         index1 = ia(node1) + iwork(node1)
    2720:         iwork(node1) = iwork(node1) + 1
    2721:         ja(index1) = node2
    2722:         index2 = ia(node2) + iwork(node2)
    2723:         iwork(node2) = iwork(node2) + 1
    2724:         ja(index2) = node1
    2725: #else
    2726:         if (node1 .le. nnodes) then
    2727:           index1 = ia(node1) + iwork(node1)
    2728:           iwork(node1) = iwork(node1) + 1
    2729:           ja(index1) = node2
    2730:         endif
    2731:         if (node2 .le. nnodes) then
    2732:           index2 = ia(node2) + iwork(node2)
    2733:           iwork(node2) = iwork(node2) + 1
    2734:           ja(index2) = node1
    2735:         endif
    2736: #endif
    2737:  1030 continue
    2738: !
    2739: ! Now lets sort all our "bins" and get the correct one on the diagonal
    2740: !
    2741:       do 1050 i = 1,nnodes
    2742:         istart = ia(i)
    2743:         iend   = ia(i+1) - 1
    2744: !       write(9,200)i,istart,iend
    2745: ! 200   format(1h ,'Sorting ',i6,' istart iend = ',i6,1x,i6)
    2746:         call SORTER(istart,iend,ja,i)
    2747:  1050 continue
    2748: !
    2749: ! Now get the "fhelp" array which will assist in assembling
    2750: ! the flux Jacobians into the correct location in the alu array
    2751: !
    2752: !     write(90,*) 'fhelp array'
    2753: !     do 1060 i = 1,nedge
    2754: !       node1 = eptr(i,1)
    2755: !       node2 = eptr(i,2)
    2756: !
    2757: ! First take care of node1
    2758: !
    2759: !       idiag = iau(node1)
    2760: !
    2761: ! If the offdiagonal term is ordered later in the ja array
    2762: !
    2763: !       if (node2.gt.node1)then
    2764: !        jstart = idiag + 1
    2765: !        jend   = ia(node1+1) - 1
    2766: !       else
    2767: !        jstart = ia(node1)
    2768: !        jend   = idiag -1
    2769: !       end if
    2770: !
    2771: !        do 1070 j = jstart,jend
    2772: !          if (ja(j).eq.node2)fhelp(i,1) = j
    2773: !1070    continue
    2774: !
    2775: !
    2776: ! Now take care of node2
    2777: !
    2778: !       idiag = iau(node2)
    2779: !
    2780: ! If the offdiagonal term is ordered later in the ja array
    2781: !
    2782: !       if (node1.gt.node2)then
    2783: !        jstart = idiag + 1
    2784: !        jend   = ia(node2+1) - 1
    2785: !       else
    2786: !        jstart = ia(node2)
    2787: !        jend   = idiag -1
    2788: !       end if
    2789: !
    2790: !       do 1080 j = jstart,jend
    2791: !         if (ja(j).eq.node1)fhelp(i,2) = j
    2792: !1080   continue
    2793: !       write(90,*) i,fhelp(i,1),fhelp(i,2)
    2794: !1060 continue
    2795: !      close(90)
    2796: !
    2797:       return
    2798:       end
    2799: !
    2800: !
    2801: !===================================================================
    2802: !
    2803: ! Sort each of our bins
    2804: !
    2805: !===================================================================
    2806:       subroutine SORTER(istart,iend,ja,inode)
    2807:       implicit none
    2808: #include <petsc/finclude/petscsys.h>
    2809:       integer istart,iend,inode
    2810:       integer ja(1)
    2811: !
    2812:       integer min,minsave,jsave,i,j
    2813:       do 1000 i = istart,iend
    2814:         min = ja(i)
    2815:         minsave = ja(i)
    2816:         jsave = i
    2817:         do 1010 j = i+1,iend
    2818:           if (ja(j).lt.min)then
    2819:             min = ja(j)
    2820:             jsave = j
    2821:           end if
    2822:  1010   continue
    2823:         ja(i) = min
    2824:         ja(jsave) = minsave
    2825: !       if (ja(i).eq.inode)iau(inode) = i
    2826:  1000 continue
    2827: !
    2828:       return
    2829:       end