Actual source code: ex74f.F90

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1:       program radhyd
  2: ! "$Id: ex4f.F,v 1.39 1999/03/10 19:29:25 Vince Mousseau $";
  3: !
  4: !  Description: This example solves a nonlinear system on 1 processor with SNES.
  5: !  We solve the Euler equations in one dimension.
  6: !  The command line options include:
  7: !    -dt <dt>, where <dt> indicates time step
  8: !    -mx <xg>, where <xg> = number of grid points in the x-direction
  9: !    -nstep <nstep>, where <nstep> = number of time steps
 10: !    -debug <ndb>, where <ndb> = 0) no debug 1) debug
 11: !    -pcnew <npc>, where <npc> = 0) no preconditioning 1) rad preconditioning
 12: !    -probnum <probnum>, where <probnum> = 1) cyclic Riesner 2) dam break
 13: !    -ihod <ihod>, where <ihod> = 1) upwind 2) quick 3) godunov
 14: !    -ientro <ientro>, where <ientro> = 0) basic 1) entropy fix 2) hlle
 15: !    -theta <theta>, where <theta> = 0-1 0-explicit 1-implicit
 16: !    -hnaught <hnaught>, where <hnaught> = height of left side
 17: !    -hlo <hlo>, where <hlo> = hieght of right side
 18: !    -ngraph <ngraph>, where <ngraph> = number of time steps between graphics
 19: !    -damfac <damfac>, where <damfac> = fractional downward change in hight
 20: !    -dampit <ndamp>, where <ndamp> = 1 turn damping on
 21: !    -gorder <gorder>, where <gorder> = spatial oerder of godunov
 22: !
 23: !
 24: !  --------------------------------------------------------------------------
 25: !
 26: ! Shock tube example
 27: !
 28: !  In this example the application context is a Fortran integer array:
 29: !      ctx(1) = shell preconditioner pressure matrix contex
 30: !      ctx(2) = semi implicit pressure matrix
 31: !      ctx(4) = xold  - old time values need for time advancement
 32: !      ctx(5) = mx    - number of control volumes
 33: !      ctx(6) = N     - total number of unknowns
 34: !
 35: !  --------------------------------------------------------------------------

 37:       implicit none

 39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 40: !                    Include files
 41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42: !
 43:  #include finclude/petscdef.h
 44:  #include finclude/petsc.h

 46: #include "ex74fcomd.h"
 47: #include "ex74ftube.h"

 49: !
 50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51: !                   Variable declarations
 52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53: !
 54: !  Variables:
 55: !     snes        - nonlinear solver
 56: !     x,r         - solution, residual vectors
 57: !     J           - Jacobian matrix
 58: !     its         - iterations for convergence
 59: !     dt          - time step size
 60: !     draw        - drawing context
 61: !
 62:       PetscFortranAddr   ctx(6)
 63:       integer            rank,size,nx,ny
 64:       SNES               snes
 65:       KSP                ksp
 66:       PC                 pc
 67:       Vec                x,r
 68:       PetscViewer        view0, view1, view2, view3, view4
 69:       Mat                Psemi
 70:       integer            flg, N, ierr, ngraph
 71:       integer            nstep, ndt,  i
 72:       integer            its, lits, totits, totlits
 73:       integer            ndb, npc, ndamp, nwilson, ndtcon
 74:       double precision   plotim

 76:       double precision krtol,katol,kdtol
 77:       double precision natol,nrtol,nstol
 78:       integer  kmit,nmit,nmf


 81: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 82: !  MUST be declared as external.

 84:       external FormFunction, FormInitialGuess,FormDt,PCRadSetUp, PCRadApply, FormGraph,FormDampit
 85:       double precision eos

 87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88: !  Initialize program
 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 91:       open (unit=87,file='Dt.out',status='unknown')


 94: !  start PETSc
 95: !
 96:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 97:       call PetscOptionsSetValue('-snes_mf','true',ierr)
 98:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 99:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

101:       if (size .ne. 1) then
102:          if (rank .eq. 0) then
103:             write(6,*) 'This is a uniprocessor example only!'
104:          endif
105:          SETERRQ(PETSC_COMM_WORLD,1,0,' ')
106:       endif

108: !  Initialize problem parameters

110:       debug       = .false.
111:       dampit      = .false.
112:       wilson      = .true.
113:       dtcon       = .true.
114:       pcnew       = .true.
115:       dtmax       = 1.0d+2
116:       dtmin       = 1.00d-12
117:       dt          = 1.0d-2
118:       mx          = 100
119:       nstep       = 50
120:       probnum     = 1
121:       gorder      = 1

123:       tfinal      = 1.0d+0
124:       tplot       = 0.2d+0
125:       dtgrow      = 1.05d+0
126:       tcscal      = 0.5d+0
127:       hcscal      = 0.5d+0

129:       ihod = 3
130:       ientro = 1
131:       theta = 1.00d+0
132:       pi = 3.14159d+0

134:       zero = 0.0
135:       ngraph = 10

137:       ndb = 0
138:       npc = 1

140:       damfac = 0.9d+0

142:       gamma = 1.25d+0
143:       csubv = 1.0d+0 / (gamma - 1.0d+0)

145:       v1 = 0.0d+0
146:       v4 = 0.0d+0

148:       e1 = 1.0d+0
149:       e4 = 1.0d+0

151:       r1 = 1.0d+0
152:       r4 = 2.0d+0

154:       ru1 = r1 * v1
155:       ru4 = r4 * v4

157:       et1 = r1 * ( (0.5d+0 * v1 * v1) + e1 )
158:       et4 = r4 * ( (0.5d+0 * v4 * v4) + e4 )

160:       p1 = eos(r1,ru1,et1)
161:       p4 = eos(r4,ru4,et4)

163:       a1 = sqrt(gamma*p1/r1)
164:       a4 = sqrt(gamma*p4/r4)

166:       erg0   = 1.0d+2
167:       kappa0 = 1.0d+0
168:       kappaa = -2.0d+0
169:       kappab = 13.0d+0 / 2.0d+0

171: !
172: !  load the command line options
173: !
174:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
175:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
176:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nstep',nstep,flg,ierr)
177:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-debug',ndb,flg,ierr)
178:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-pcnew',npc,flg,ierr)
179:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ihod',ihod,flg,ierr)
180:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ientro',ientro,flg,ierr)
181:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',theta,flg,ierr)
182:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ngraph',ngraph,flg,ierr)
183:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-damfac',damfac,flg,ierr)
184:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dampit',ndamp,flg,ierr)
185:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER, '-wilson',nwilson,flg,ierr)
186:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-gorder',gorder,flg,ierr)
187:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER, '-probnum',probnum,flg,ierr)
188:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-kappa0',kappa0,flg,ierr)
189:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-erg0',erg0,flg,ierr)
190:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dtcon',ndtcon,flg,ierr)
191:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tfinal',tfinal,flg,ierr)
192:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tplot',tplot,flg,ierr)
193:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dtgrow',dtgrow,flg,ierr)
194:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tcscal',tcscal,flg,ierr)
195:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-hcscal',hcscal,flg,ierr)
196:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-dtmax',dtmax,flg,ierr)

198:       if (ndamp .eq. 1) then
199:          dampit = .true.
200:       endif

202:       if (nwilson .eq. 0) then
203:          wilson = .false.
204:       endif

206:       if (ndb .eq. 1) then
207:          debug = .true.
208:       endif

210:       if (npc .eq. 0) then
211:          pcnew = .false.
212:       endif

214:       if (ndtcon .eq. 0) then
215:          dtcon = .false.
216:       endif

218: !CVAM  if (dt .ge. dtmax .or. dt .le. dtmin) then
219: !CVAM     if (rank .eq. 0) write(6,*) 'DT is out of range'
220: !CVAM     SETERRA(1,0,' ')
221: !CVAM  endif

223:       N       = mx*neq

225:       ctx(5) = mx
226:       ctx(6) = N

228:       if (debug) then
229:         write(*,*) 'mx = ',mx
230:       endif



234: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: !  Create nonlinear solver context
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

238:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

240: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: !  Create vector data structures; set function evaluation routine
242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

244:       call MatCreateSeqAIJ(PETSC_COMM_WORLD,mx,mx,10,PETSC_NULL_INTEGER,ctx(2),ierr)

246:       if (debug) then
247:         call MatGetSize(ctx(2),nx,ny,ierr)
248:         write(*,*) 'number of rows = ',nx,' number of col = ',ny
249:       endif
250: !
251: !  full size vectors
252: !
253:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,x,ierr)
254:       call VecSetFromOptions(x,ierr)
255:       call VecDuplicate(x,r,ierr)
256:       call VecDuplicate(x,ctx(4),ierr)
257: !
258: ! set grid
259: !
260:       dx = 2.0d+0/dfloat(mx)
261:       xl0 = -1.0d+0 -(0.5d+0 * dx)

263:       if (debug) then
264:         write(*,*) 'dx = ',dx
265:       endif
266: 

268: !  Set function evaluation routine and vector.  Whenever the nonlinear
269: !  solver needs to evaluate the nonlinear function, it will call this
270: !  routine.
271: !   - Note that the final routine argument is the user-defined
272: !     context that provides application-specific data for the
273: !     function evaluation routine.

275:       call SNESSetFunction(snes,r,FormFunction,ctx,ierr)

277: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: !  Customize nonlinear solver; set runtime options
279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

281: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)

283:       call SNESSetFromOptions(snes,ierr)
284: !
285: !  set the line search function to damp the newton update.
286: !
287:       if (dampit) then
288:         call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
289:       endif
290: !
291: ! get the linear solver info from the nonlinear solver
292: !

294:       call SNESGetKSP(snes,ksp,ierr)
295:       call KSPGetPC(ksp,pc,ierr)

297:       call KSPGetTolerances(ksp,krtol,katol,kdtol,kmit,ierr)
298:       call SNESGetTolerances(snes,natol,nrtol,nstol,nmit,nmf,ierr)

300:       write(*,*)
301:       write(*,*)
302:       write(*,*) 'Linear solver'
303:       write(*,*)
304:       write(*,*) 'rtol = ',krtol
305:       write(*,*) 'atol = ',katol
306:       write(*,*) 'dtol = ',kdtol
307:       write(*,*) 'maxits = ',kmit
308:       write(*,*)
309:       write(*,*)
310:       write(*,*) 'Nonlinear solver'
311:       write(*,*)
312:       write(*,*) 'rtol = ',nrtol
313:       write(*,*) 'atol = ',natol
314:       write(*,*) 'stol = ',nstol
315:       write(*,*) 'maxits = ',nmit
316:       write(*,*) 'max func = ',nmf
317:       write(*,*)
318:       write(*,*)

320: !
321: !  Build shell based preconditioner if flag set
322: !
323:       if (pcnew) then
324:         call PCSetType(pc,PCSHELL,ierr)
325:         call PCShellSetContext(pc,ctx,ierr)
326:         call PCShellSetSetUpCtx(pc,PCRadSetUp,ierr)
327:         call PCShellSetApplyCtx(pc,PCRadApply,ctx,ierr)
328:       endif

330:       call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)

332: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: !  Evaluate initial guess; then solve nonlinear system.
334: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335: !
336: !  initial counters
337: !
338:       time = 0.0d+0
339:       plotim = 0.0d+0
340:       totits = 0
341:       totlits = 0

343: !  Note: The user should initialize the vector, x, with the initial guess
344: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
345: !  to employ an initial guess of zero, the user should explicitly set
346: !  this vector to zero by calling VecSet().

348:       call FormInitialGuess(x,ierr)
349: !
350: !  open a window to plot results
351: !
352:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'density',0,0,300,300,view0,ierr)
353:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'velocity',320,0,300,300,view1,ierr)
354:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'total energy',640,0,300,300,view2,ierr)
355:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'temperature',0,360,300,300,view3,ierr)
356:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'pressure',320,360,300,300,view4,ierr)
357: !
358: !  graph initial conditions
359: !
360:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)
361: !
362: !  copy x into xold
363: !
364:       call VecCopy(x,ctx(4),ierr)
365:       call FormDt(snes,x,ctx,ierr)
366: !
367: !################################
368: !
369: !  TIME STEP LOOP BEGIN
370: !
371: !################################
372: !
373:       ndt = 0

375:    10 if ( (ndt .le. nstep) .and. ((time + 1.0d-10) .lt. tfinal) ) then

377:         if (debug) then
378:           write(*,*)
379:           write(*,*) 'start of time loop'
380:           write(*,*)
381:           write(*,*) 'ndt = ',ndt
382:           write(*,*) 'nstep = ',nstep
383:           write(*,*) 'time = ',time
384:           write(*,*) 'tfinal = ',tfinal
385:           write(*,*)
386:         endif

388:         ndt = ndt + 1
389: !
390: !  increment time
391: !
392:         time = time + dt
393:         plotim = plotim + dt
394: !
395: !  call the nonlinear solver
396: !
397:         call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
398: !
399: !  get the number of linear iterations used by the nonlinear solver
400: !
401:         call SNESGetLinearSolveIterations(snes,lits,ierr)
402:         call SNESGetIterationNumber(snes,its,ierr)

404:         if (debug) then
405:            write(*,*) 'in radhyd ',ndt,'x'
406:            call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
407:         endif
408: !
409: !  increment the counters
410: !
411:         totits = totits + its
412:         totlits = totlits + lits
413: !
414: !  Compute new time step
415: !
416:           call FormDt(snes,x,ctx,ierr)
417: !
418: !  Draw contour plot of solution
419: !
420:         if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
421: 
422:            plotim = plotim - tplot


425:         if (rank .eq. 0) then
426:            write(6,100) totits,totlits,ndt,dt,time
427:         endif
428:   100   format('Newt = ',i7,' lin =',i7,' step =',i7,' dt = ',e8.3,' time = ',e10.4)
429: !
430: !  graph state conditions
431: !
432:           call FormGraph(x,view0,view1,view2,view3,view4,ierr)

434:         endif
435: !
436: ! copy x into xold
437: !
438:         call VecCopy(x,ctx(4),ierr)


441:         goto 10

443:       endif
444: !
445: !################################
446: !
447: !  TIME STEP LOOP END
448: !
449: !################################
450: !

452: !
453: !  graph final conditions
454: !
455:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)


458:       write(*,*)
459:       write(*,*)
460:       write(*,*) 'total Newton iterations = ',totits
461:       write(*,*) 'total linear iterations = ',totlits
462:       write(*,*) 'Average Newton per time step = ', dble(totits)/dble(ndt)
463:       write(*,*) 'Average Krylov per Newton = ', dble(totlits)/dble(totits)
464:       write(*,*)
465:       write(*,*)

467: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
468: !  Free work space.  All PETSc objects should be destroyed when they
469: !  are no longer needed.
470: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


473:       call MatDestroy(ctx(2),ierr)
474:       call VecDestroy(x,ierr)
475:       call VecDestroy(ctx(4),ierr)
476:       call VecDestroy(r,ierr)
477:       call SNESDestroy(snes,ierr)
478:       call PetscViewerDestroy(view0,ierr)
479:       call PetscViewerDestroy(view1,ierr)
480:       call PetscViewerDestroy(view2,ierr)
481:       call PetscViewerDestroy(view3,ierr)
482:       call PetscViewerDestroy(view4,ierr)

484:       call PCDestroy(ctx(1),ierr)

486:       call PetscFinalize(ierr)

488:       close(87)

490:       stop
491:       end
492:       subroutine ApplicationDampit(x,deltx,w,ierr)
493: ! ---------------------------------------------------------------------
494: !
495: !  ApplicationDampit - Damps the newton update, called by
496: !  the higher level routine FormDampit().
497: !
498: !  Input Parameter:
499: !  x    - current iterate
500: !  deltx   - update
501: !
502: !  Output Parameters:
503: !  x    - new iterate
504: !  ierr - error code
505: !
506: !  Notes:
507: !  This routine only damps the density.  May want to add energy
508: !  in the future
509: !

511:       implicit none

513: !  Common blocks:
514: #include "ex74fcomd.h"

516: !  Input/output variables:
517:       PetscScalar   x(mx*neq), deltx(mx*neq), w(mx*neq)
518:       integer  ierr

520: !  Local variables:
521:       double precision facmin, fac, newx, xmin, se, dse
522:       double precision u,en,rn,run
523:       integer  i, jr, jru, je

525:       0

527:       if (debug) then
528:         write(*,*) 'begin damping'
529:         do i = 1,mx*neq
530:           write(*,*)'x(',i,') = ',x(i)
531:         enddo
532:         write(*,*)
533:         do i = 1,mx*neq
534:           write(*,*)'deltx(',i,') = ',deltx(i)
535:         enddo
536:       endif

538:       facmin = 1.0d+0
539: !
540: !  set the scale factor
541: !
542:       do i=1,mx
543: !
544: !  set pointers
545: !
546:         jr  = (neq*i) - 2
547:         jru = (neq*i) - 1
548:         je  = (neq*i)
549: !
550: !  make sure dencity stayes positive
551: !
552:         newx = x(jr) - deltx(jr)
553:         xmin = damfac * x(jr)

555:         if (newx .lt. xmin) then
556:           fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
557:           if (fac .lt. facmin) then
558:             if (debug) then
559:               write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
560:             endif
561:             facmin = fac
562:           endif
563:         endif
564: !
565: !  make sure Total energy stayes positive
566: !
567:         newx = x(je) - deltx(je)
568:         xmin = damfac * x(je)

570:         if (newx .lt. xmin) then
571:           fac = (1.0d+0 - damfac)*x(je)/deltx(je)
572:           if (fac .lt. facmin) then
573:             if (debug) then
574:               write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
575:             endif
576:             facmin = fac
577:           endif
578:         endif
579: !
580: !  make sure specific internal  energy stayes positive
581: !
582: 
583:         u = x(jru)/x(jr)
584:         se = (x(je)/x(jr)) - (0.5d+0 * u * u)

586:         en = x(je) - deltx(je)
587:         rn = x(jr) - deltx(jr)
588:         run = x(jru) - deltx(jru)

590:         dse = se - ( (en/rn) - (0.5d+0 * (run/rn) * (run/rn)) )


593:         newx = se - dse
594:         xmin = damfac * se

596:         if (newx .lt. xmin) then
597:           fac = (1.0d+0 - damfac) * se / dse
598:           if (fac .lt. facmin) then
599:             if (debug) then
600:               write(*,*) 'se',i, damfac,facmin,fac,se,dse
601:             endif
602:             facmin = fac
603:           endif
604:         endif

606:       enddo
607: !
608: ! write out warning
609: !
610:       if (facmin .lt. 1.0d+0) then
611:         write(*,*) 'facmin = ',facmin, damfac,time
612: !
613: !  scale the vector
614: !
615:         do i=1,neq*mx
616:           w(i) = x(i) - (facmin * deltx(i))
617:         enddo
618:       else
619:         do i=1,neq*mx
620:           w(i) = x(i) -  deltx(i)
621:         enddo
622:       endif

624:       if (debug) then
625:         write(*,*) 'end damping'
626:         do i = 1,mx*neq
627:            write(*,*) 'w(',i,') = ',w(i)
628:         enddo
629:       endif

631:       return
632:       end
633:       subroutine ApplicationDt(x,xold,ierr)
634: ! ---------------------------------------------------------------------
635: !
636: !  ApplicationDt - compute CFL numbers. Called by
637: !  the higher level routine FormDt().
638: !
639: !  Input Parameter:
640: !  x    - local vector data
641: !
642: !  Output Parameters:
643: !  ierr - error code
644: !
645: !  Notes:
646: !  This routine uses standard Fortran-style computations over a 2-dim array.
647: !

649:       implicit none

651: !  Common blocks:
652: #include "ex74fcomd.h"
653: #include "ex74ftube.h"

655: !  Input/output variables:
656:       PetscScalar   x(mx*neq), xold(mx*neq)
657:       integer  ierr

659: !  Local variables:
660:       integer  i, jr, jru, je
661: !
662: ! new
663: !
664:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw
665:       double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
666:       double precision tcflg, hcflg, dtt, dth
667:       double precision te, tp, tw
668:       double precision ue, up, uw
669:       double precision see, sep, sew
670: !
671: ! old
672: !
673:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
674:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
675:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
676:       double precision veloe,  velop,  velow
677:       double precision uop, seop, top
678:       double precision dtold, dttype
679: !
680: !  functions
681: !
682:       double precision eos

684:       dtold = dt

686:       0

688:       if (debug) then
689:         write(*,*) 'in start dt'
690:         do i = 1,mx*neq
691:           write(*,*)'x(',i,') = ',x(i)
692:         enddo
693:         write(*,*) 'tfinal = ',tfinal
694:         write(*,*) 'time = ',time
695:         write(*,*) 'dt = ',dt
696:         write(*,*) 'dtmax = ',dtmax
697:       endif

699:       sndp = -1.0d+20
700:       vradn = 0.0d+0
701:       vradd = 0.0d+0

703: !
704: !################################
705: !
706: !  loop over all cells begin
707: !
708: !################################
709: !
710:       do i=1,mx
711: !
712: !  set pointers
713: !
714:         jr  = (neq*i) - 2
715:         jru = (neq*i) - 1
716:         je  = (neq*i)
717: !
718: !
719: !  set scalars
720: !
721:         call Setpbc(i,x, rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww,vele,  velp,  velw)
722: !
723: ! compute temperatures
724: !
725:         uw = rhouw / rhow
726:         up = rhoup / rhop
727:         ue = rhoue / rhoe

729:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
730:         sep = (ergp/rhop) - (0.5d+0 * up * up)
731:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

733:         te  = see / csubv
734:         tp  = sep / csubv
735:         tw  = sew / csubv
736: !
737: ! compute old temperature
738: !

740:         call Setpbc(i,xold,rhooee,  rhooe,  rhoop,  rhoow,  rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow)

742:         call Setpbcn(rhooee,  rhooe,  rhoop,  rhoow,  rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww,veloe,  velop,  velow, i)

744:         uop = rhouop / rhoop

746:         seop = (ergop/rhoop) - (0.5d+0 * uop * uop)

748:         top  = seop / csubv
749: !
750: !  compute thermal cfl
751: !
752:         vradn = vradn + (abs(tp - top)/dt)
753:         vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
754: !
755: !  compute hydro cfl
756: !

758:         pressp  = eos(rhop, rhoup, ergp)
759:         sndp = max(sndp,sqrt( (gamma*pressp) / rhop ))

761:       enddo
762: !
763: !################################
764: !
765: !  loop over all cells end
766: !
767: !################################
768: !

770:       vrad = vradn / vradd

772:       tcfl = (vrad * dt) / dx
773:       hcfl = (sndp * dt) / dx

775:       dtt = max(dx/vrad,1.0d-7)
776:       dtt = tcscal * dtt

778:       dth = hcscal * dx / sndp

780:       if (.not. dtcon) then
781:         dt = min (dth,dtt,dt*dtgrow)
782:         if (dt .lt. dtmin) then
783:            dt = dtmin
784:         endif
785:         if (dt .gt. dtmax) then
786:            dt = dtmax
787:         endif
788:         if ( (time + dt) .gt. tfinal) then
789:            dt = tfinal - time
790:         endif

792:         if (dt .eq. dth) then
793:            dttype = 1.0d+0
794:         elseif (dt .eq. dtt) then
795:            dttype = 2.0d+0
796:         elseif (dt .eq. dtold*dtgrow) then
797:            dttype = 3.0d+0
798:         elseif (dt .eq. dtmax) then
799:            dttype = 4.0d+0
800:         elseif (dt .eq. dtmin) then
801:            dttype = 5.0d+0
802:         elseif (dt .eq. tfinal - time) then
803:            dttype = 6.0
804:         else
805:            dttype = -1.0d+0
806:         endif

808:       endif
809: 
810: 
811:       write(87,1000) time,dt,dth/hcscal,dtt/tcscal
812:       write(88,1000) time,dttype

814:  1000 format(4(2x,e18.9))

816:       if (debug) then
817:         write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
818:         write(*,*) 'dtt = ',dtt,' dth = ',dth
819:         write(*,*) 'tfinal = ',tfinal
820:         write(*,*) 'time = ',time
821:         write(*,*) 'dt = ',dt
822:         write(*,*) 'dtmax = ',dtmax
823:         write(*,*)
824:       endif

826:       return
827:       end
828:       subroutine ApplicationExact(x,ierr)
829: ! ---------------------------------------------------------------------
830: !
831: !  ApplicationExact - Computes exact solution, called by
832: !  the higher level routine FormExact().
833: !
834: !  Input Parameter:
835: !  x - local vector data
836: !
837: !  Output Parameters:
838: !  x -    initial conditions
839: !  ierr - error code
840: !
841: !  Notes:
842: !  This routine uses standard Fortran-style computations over a 1-dim array.
843: !

845:       implicit none

847: !  Common blocks:

849: #include "ex74fcomd.h"

851: !  Input/output variables:
852:       PetscScalar  x(mx)
853:       integer ierr

855: !  Local variables:
856:       integer  i
857:       double precision xloc
858:       PetscScalar rexact


861: !  Set parameters

863:       0

865:       do i = 1,mx

867:         xloc = xl0 + (dble(i) * dx)
868:         x(i) = rexact(xloc,time)

870:       enddo

872:       return
873:       end
874:       subroutine ApplicationFunction(x,f,xold,ierr)
875: ! ---------------------------------------------------------------------
876: !
877: !  ApplicationFunction - Computes nonlinear function, called by
878: !  the higher level routine FormFunction().
879: !
880: !  Input Parameter:
881: !  x    - local vector data
882: !
883: !  Output Parameters:
884: !  f    - local vector data, f(x)
885: !  ierr - error code
886: !
887: !  Notes:
888: !  This routine uses standard Fortran-style computations over a 2-dim array.
889: !

891:       implicit none

893: !  Common blocks:
894: #include "ex74fcomd.h"

896: !  Input/output variables:
897:       PetscScalar   x(mx*neq), f(mx*neq), xold(mx*neq)
898:       integer  ierr

900: !  Local variables:
901:       integer  i, jr, jru, je
902:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw

904:       double precision cont, energy, mom

906:       0

908:       if (debug) then
909:         write(*,*) 'in function'
910:         do i = 1,mx*neq
911:           write(*,*)'x(',i,') = ',x(i)
912:         enddo
913:       endif
914: !
915: !################################
916: !
917: !  loop over all cells begin
918: !
919: !################################
920: !
921:       do i=1,mx
922: !
923: !  set pointers
924: !
925:       jr  = (neq*i) - 2
926:       jru = (neq*i) - 1
927:       je  = (neq*i)
928: !
929: !
930: !  set scalars
931: !
932:       call Setpbc(i,x,rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw)
933: !
934: !  compute functions
935: !

937:        f(jr) = cont(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, i,xold)


940:        f(jru) = mom(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, i,xold)


943:        f(je) = energy(rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, i,xold)

945:        if (debug) then
946:          write(*,*)
947:          write(*,*) i,jr,jru,je,'res,r,ru,e'
948:          write(*,*) f(jr),f(jru),f(je)
949:          write(*,*)
950:        endif

952:       enddo
953: !
954: !################################
955: !
956: !  loop over all cells end
957: !
958: !################################
959: !

961:       if (debug) then
962:         write(*,*) 'in function'
963:         do i = 1,mx*neq
964:            write(*,*) 'f(',i,') = ',f(i)
965:         enddo
966:       endif

968:       return
969:       end
970:       subroutine ApplicationInitialGuess(x,ierr)
971: ! ---------------------------------------------------------------------
972: !
973: !  ApplicationInitialGuess - Computes initial approximation, called by
974: !  the higher level routine FormInitialGuess().
975: !
976: !  Input Parameter:
977: !  x - local vector data
978: !
979: !  Output Parameters:
980: !  x -    initial conditions
981: !  ierr - error code
982: !
983: !  Notes:
984: !  This routine uses standard Fortran-style computations over a 1-dim array.
985: !

987:       implicit none

989: !  Common blocks:

991: #include "ex74fcomd.h"
992: #include "ex74ftube.h"

994: !  Input/output variables:
995:       PetscScalar  x(mx*neq)
996:       integer ierr

998: !  Local variables:
999:       integer  i, j, jr, jru, je
1000:       double precision xloc, re, ee, ve
1001:       double precision wid, efloor
1002:       PetscScalar uexact, rexact, eexact


1005: !VAM  efloor = max(1.0d-1,1.0d-3 * erg0)
1006:       efloor = 1.0d-1
1007: !VAM  wid = max(1.0d-2,dx)
1008:       wid = 1.0d-2

1010: !  Set parameters

1012:       0

1014:       do i = 1,mx

1016:         jr  = (neq*i) - 2
1017:         jru = (neq*i) - 1
1018:         je  = (neq*i)

1020:         xloc = xl0 + (dble(i) * dx)

1022:         if (probnum .eq. 1) then
1023:            re = rexact(xloc,zero)
1024:            ve = uexact(xloc,zero)
1025:            ee = eexact(xloc,zero)
1026:         else
1027:            re = 1.0d+0
1028:            ve = 0.0d+0
1029:            ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1030:         endif

1032:         x(jr)  = re
1033:         x(jru) = re * ve
1034:         x(je)  = re * ( (0.5d+0 * ve * ve) + ee )

1036:         if (debug) then
1037:            write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1038:  100       format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1039:         endif

1041:       enddo

1043:       call exact0
1044:       call eval2
1045:       call rval2
1046:       call wval
1047:       call uval2
1048:       v3 = v2
1049:       call val3

1051:       a1 = sqrt(gamma*p1/r1)
1052:       a2 = sqrt(gamma*p2/r2)
1053:       a3 = sqrt(gamma*p3/r3)
1054:       a4 = sqrt(gamma*p4/r4)

1056:       write(*,1000) r1,r2,r3,r4
1057:       write(*,2000) p1,p2,p3,p4
1058:       write(*,3000) e1,e2,e3,e4
1059:       write(*,4000) a1,a2,a3,a4
1060:       write(*,*)

1062:  1000 format ('rhos      ',4(f12.6))
1063:  2000 format ('pressures ',4(f12.6))
1064:  3000 format ('energies  ',4(f12.6))
1065:  4000 format ('sound     ',4(f12.6))


1068:       return
1069:       end
1070:       subroutine ApplicationXmgr(x,ivar,ierr)
1071: ! ---------------------------------------------------------------------
1072: !
1073: !  ApplicationXmgr - Sets the Xmgr output called from
1074: !  the higher level routine FormXmgr().
1075: !
1076: !  Input Parameter:
1077: !  x - local vector data
1078: !
1079: !  Output Parameters:
1080: !  x -    initial conditions
1081: !  ierr - error code
1082: !
1083: !  Notes:
1084: !  This routine uses standard Fortran-style computations over a 1-dim array.
1085: !

1087:       implicit none

1089: !  Common blocks:

1091: #include "ex74fcomd.h"

1093: !  Input/output variables:
1094:       PetscScalar  x(mx)
1095:       integer ivar,ierr

1097: !  Local variables:
1098:       integer  i
1099:       double precision xloc, sum
1100:       PetscScalar rexact
1101:       integer iplotnum(5)
1102:       save iplotnum
1103:       character*8 grfile


1106:       data iplotnum / -1,-1,-1,-1,-1 /



1110: !  Set parameters

1112:       iplotnum(ivar) = iplotnum(ivar) + 1
1113:       0

1115:       if (ivar .eq. 1) then
1116:          write(grfile,4000) iplotnum(ivar)
1117:  4000    format('Xmgrr',i3.3)
1118:       elseif (ivar .eq. 2) then
1119:          write(grfile,5000) iplotnum(ivar)
1120:  5000    format('Xmgru',i3.3)
1121:       elseif (ivar .eq. 3) then
1122:          write(grfile,6000) iplotnum(ivar)
1123:  6000    format('Xmgre',i3.3)
1124:       elseif (ivar .eq. 4) then
1125:          write(grfile,7000) iplotnum(ivar)
1126:  7000    format('Xmgrt',i3.3)
1127:       else
1128:          write(grfile,8000) iplotnum(ivar)
1129:  8000    format('Xmgrp',i3.3)
1130:       endif

1132:       open (unit=44,file=grfile,status='unknown')

1134:       do i = 1,mx

1136:         xloc = xl0 + (dble(i) * dx)
1137:         if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1138:           write(44,1000) xloc, x(i), rexact(xloc,time)
1139:         else
1140:           write(44,1000) xloc, x(i)
1141:         endif

1143:       enddo

1145:  1000 format(3(e18.12,2x))
1146:       close(44)

1148:       if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1149:         sum = 0.0d+0
1150:         do i = 1,mx
1151:            xloc = xl0 + (dble(i) * dx)
1152:            sum = sum + (x(i) - rexact(xloc,time)) ** 2
1153:         enddo
1154:         sum = sqrt(sum)

1156:         write(*,*)
1157:         write(*,*)  'l2norm of the density error is',sum
1158:         write(*,*)
1159:       endif


1162:       return
1163:       end
1164:       subroutine FormDampit(snes,ctx,x,f,g,y,w, fnorm,ynorm,gnorm,flag,ierr)
1165: ! ---------------------------------------------------------------------
1166: !
1167: !  FormDampit - damps the Newton update
1168: !
1169: !  Input Parameters:
1170: !  snes  - the SNES context
1171: !  x     - current iterate
1172: !  f     - residual evaluated at x
1173: !  y     - search direction (containes new iterate on output)
1174: !  w     - work vector
1175: !  fnorm - 2-norm of f
1176: !
1177: !  In this example the application context is a Fortran integer array:
1178: !      ctx(1) = shell preconditioner pressure matrix contex
1179: !      ctx(2) = semi implicit pressure matrix
1180: !      ctx(4) = xold  - old time values need for time advancement
1181: !      ctx(5) = mx    - number of control volumes
1182: !      ctx(6) = N     - total number of unknowns
1183: !
1184: !  Output Parameter:
1185: !  g     - residual evaluated at new iterate y
1186: !  y     - new iterate (contains search direction on input
1187: !  gnorm - 2-norm of g
1188: !  ynorm - 2-norm of search length
1189: !  flag  - set to 0 if the line search succeeds; -1 on failure
1190: !
1191: !  Notes:
1192: !  This routine serves as a wrapper for the lower-level routine
1193: !  "ApplicationDampit", where the actual computations are
1194: !  done using the standard Fortran style of treating the local
1195: !  vector data as a multidimensional array over the local mesh.
1196: !  This routine merely accesses the local vector data via
1197: !  VecGetArray() and VecRestoreArray().
1198: !
1199:       implicit none

1201:  #include finclude/petsc.h

1203: !  Input/output variables:
1204:       SNES             snes
1205:       Vec              x, f, g, y, w
1206:       PetscFortranAddr ctx(*)
1207:       PetscScalar           fnorm, ynorm, gnorm
1208:       integer          ierr, flag

1210: !  Common blocks:

1212: #include "ex74fcomd.h"

1214: !  Local variables:

1216: !  Declarations for use with local arrays:
1217:       PetscScalar      lx_v(0:1), ly_v(0:1), lw_v(0:1)
1218:       PetscOffset lx_i, ly_i , lw_i

1220: !
1221: !  set ynorm
1222: !
1223:       call VecNorm(y,NORM_2,ynorm,ierr)
1224: !
1225: !  copy x to w
1226: !
1227:       call VecCopy(x,w,ierr)
1228: !
1229: ! get pointers to x, y, w
1230: !

1232:       call VecGetArray(x,lx_v,lx_i,ierr)
1233:       call VecGetArray(y,ly_v,ly_i,ierr)
1234:       call VecGetArray(w,lw_v,lw_i,ierr)
1235: !
1236: !  Compute Damping
1237: !
1238:       call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1239: !
1240: !  Restore vectors x, y, w
1241: !
1242:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1243:       call VecRestoreArray(y,ly_v,ly_i,ierr)
1244:       call VecRestoreArray(w,lw_v,lw_i,ierr)
1245: !
1246: !  copy w to y
1247: !
1248:       call VecCopy(w,y,ierr)
1249: !
1250: !  compute new residual
1251: !
1252:       call SNESComputeFunction(snes,y,g,ierr)
1253:       call VecNorm(g,NORM_2,gnorm,ierr)
1254:       flag = 0

1256:       if (debug) then
1257:          write(*,*) 'form damp ynorm = ',ynorm
1258:          write(*,*) 'gnorm = ',gnorm
1259:       endif

1261:       return
1262:       end
1263:       subroutine FormDt(snes,x,ctx,ierr)
1264: ! ---------------------------------------------------------------------
1265: !
1266: !  FormDt - Compute CFL numbers
1267: !
1268: !  Input Parameters:
1269: !  snes  - the SNES context
1270: !  x     - input vector
1271: !
1272: !  In this example the application context is a Fortran integer array:
1273: !      ctx(1) = shell preconditioner pressure matrix contex
1274: !      ctx(2) = semi implicit pressure matrix
1275: !      ctx(4) = xold  - old time values need for time advancement
1276: !      ctx(5) = mx    - number of control volumes
1277: !      ctx(6) = N     - total number of unknowns
1278: !
1279: !
1280: !  Notes:
1281: !  This routine serves as a wrapper for the lower-level routine
1282: !  "ApplicationDt", where the actual computations are
1283: !  done using the standard Fortran style of treating the local
1284: !  vector data as a multidimensional array over the local mesh.
1285: !  This routine merely accesses the local vector data via
1286: !  VecGetArray() and VecRestoreArray().
1287: !
1288:       implicit none

1290:  #include finclude/petsc.h

1292: !  Input/output variables:
1293:       SNES             snes
1294:       Vec              x
1295:       PetscFortranAddr ctx(*)
1296:       integer          ierr

1298: !  Common blocks:

1300: #include "ex74fcomd.h"

1302: !  Local variables:

1304: !  Declarations for use with local arrays:
1305:       PetscScalar      lx_v(0:1)
1306:       PetscOffset lx_i
1307:       PetscScalar      lxold_v(0:1)
1308:       PetscOffset lxold_i

1310: !
1311: ! get pointers to x, xold
1312: !

1314:       call VecGetArray(x,lx_v,lx_i,ierr)
1315:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1316: !
1317: !  Compute function
1318: !
1319:       call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1320: !
1321: !  Restore vectors x, xold
1322: !
1323:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1324:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)

1326:       return
1327:       end
1328:       subroutine FormExact(x,ierr)
1329: ! ---------------------------------------------------------------------
1330: !
1331: !  FormExact - Forms exact solution
1332: !
1333: !  Input Parameter:
1334: !  x - vector
1335: !
1336: !  Output Parameters:
1337: !  x - vector
1338: !  ierr - error code
1339: !
1340: !  Notes:
1341: !  This routine serves as a wrapper for the lower-level routine
1342: !  "ApplicationExact", where the actual computations are
1343: !  done using the standard Fortran style of treating the local
1344: !  vector data as a multidimensional array over the local mesh.
1345: !  This routine merely accesses the local vector data via
1346: !  VecGetArray() and VecRestoreArray().
1347: !
1348:       implicit none

1350:  #include finclude/petsc.h

1352: !  Input/output variables:
1353:       Vec      x
1354:       integer  ierr

1356: !  Declarations for use with local arrays:
1357:       PetscScalar      lx_v(0:1)
1358:       PetscOffset lx_i

1360:       0

1362: !
1363: !  get a pointer to x
1364: !
1365:       call VecGetArray(x,lx_v,lx_i,ierr)
1366: !
1367: !  Compute initial guess
1368: !
1369:       call ApplicationExact(lx_v(lx_i),ierr)
1370: !
1371: !  Restore vector x
1372: !
1373:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1375:       return
1376:       end
1377:       subroutine FormFunction(snes,x,f,ctx,ierr)
1378: ! ---------------------------------------------------------------------
1379: !
1380: !  FormFunction - Evaluates nonlinear function, f(x).
1381: !
1382: !  Input Parameters:
1383: !  snes  - the SNES context
1384: !  x     - input vector
1385: !
1386: !  In this example the application context is a Fortran integer array:
1387: !      ctx(1) = shell preconditioner pressure matrix contex
1388: !      ctx(2) = semi implicit pressure matrix
1389: !      ctx(4) = xold  - old time values need for time advancement
1390: !      ctx(5) = mx    - number of control volumes
1391: !      ctx(6) = N     - total number of unknowns
1392: !
1393: !  Output Parameter:
1394: !  f     - vector with newly computed function
1395: !
1396: !  Notes:
1397: !  This routine serves as a wrapper for the lower-level routine
1398: !  "ApplicationFunction", where the actual computations are
1399: !  done using the standard Fortran style of treating the local
1400: !  vector data as a multidimensional array over the local mesh.
1401: !  This routine merely accesses the local vector data via
1402: !  VecGetArray() and VecRestoreArray().
1403: !
1404:       implicit none

1406:  #include finclude/petsc.h

1408: !  Input/output variables:
1409:       SNES             snes
1410:       Vec              x, f
1411:       PetscFortranAddr ctx(*)
1412:       integer          ierr

1414: !  Common blocks:

1416: #include "ex74fcomd.h"

1418: !  Local variables:

1420: !  Declarations for use with local arrays:
1421:       PetscScalar      lx_v(0:1), lf_v(0:1)
1422:       PetscOffset lx_i, lf_i
1423:       PetscScalar      lxold_v(0:1)
1424:       PetscOffset lxold_i

1426: !
1427: ! get pointers to x, f, xold
1428: !

1430:       call VecGetArray(x,lx_v,lx_i,ierr)
1431:       call VecGetArray(f,lf_v,lf_i,ierr)
1432:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1433: !
1434: !  Compute function
1435: !
1436:       call ApplicationFunction(lx_v(lx_i),lf_v(lf_i), lxold_v(lxold_i),ierr)
1437: !
1438: !  Restore vectors x, f, xold
1439: !
1440:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1441:       call VecRestoreArray(f,lf_v,lf_i,ierr)
1442:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1443: !
1444: ! something to do with profiling
1445: !
1446:       call PetscLogFlops(11*mx,ierr)

1448:       return
1449:       end
1450:       subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1451: ! ---------------------------------------------------------------------
1452: !
1453: !  FormGraph - Forms Graphics output
1454: !
1455: !  Input Parameter:
1456: !  x - vector
1457: !  time - scalar
1458: !
1459: !  Output Parameters:
1460: !  ierr - error code
1461: !
1462: !  Notes:
1463: !  This routine serves as a wrapper for the lower-level routine
1464: !  "ApplicationXmgr", where the actual computations are
1465: !  done using the standard Fortran style of treating the local
1466: !  vector data as a multidimensional array over the local mesh.
1467: !  This routine merely accesses the local vector data via
1468: !  VecGetArray() and VecRestoreArray().
1469: !
1470:       implicit none

1472:  #include finclude/petsc.h

1474: #include "ex74fcomd.h"
1475: #include "ex74ftube.h"

1477: !  Input/output variables:
1478:       Vec      x
1479:       integer  ierr

1481: !  Declarations for use with local arrays:
1482:       IS                 rfrom, rto, rufrom, ruto, efrom, eto
1483:       Vec                rval
1484:       Vec                uval
1485:       Vec                ruval
1486:       Vec                eval
1487:       Vec                seval
1488:       Vec                pval
1489:       Vec                kval
1490:       Vec                tval
1491:       Vec                steval
1492:       VecScatter         scatter
1493:       PetscViewer             view0, view1, view2, view3, view4
1494:       double precision minus1, l2err, gm1, csubvi


1497:       csubvi = 1.0d+0 / csubv
1498:       gm1 = gamma - 1.0d+0
1499:       0
1500:       minus1 = -1.0d+0
1501: !
1502: !  graphics vectors
1503: !
1504:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,rval,ierr)
1505:       call VecSetFromOptions(rval,ierr)
1506:       call VecDuplicate(rval,uval,ierr)
1507:       call VecDuplicate(rval,ruval,ierr)
1508:       call VecDuplicate(rval,eval,ierr)
1509:       call VecDuplicate(rval,seval,ierr)
1510:       call VecDuplicate(rval,pval,ierr)
1511:       call VecDuplicate(rval,kval,ierr)
1512:       call VecDuplicate(rval,tval,ierr)
1513:       call VecDuplicate(rval,steval,ierr)
1514: !
1515: ! create index sets for vector scatters
1516: !
1517:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1518:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  rto,   ierr)
1519:       call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1520:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  ruto,  ierr)
1521:       call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1522:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  eto,   ierr)

1524: !
1525: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1526: !
1527: !
1528: !  load rval from x
1529: !
1530:       call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1531:       call VecScatterBegin(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1532:       call VecScatterEnd(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1533:       call VecScatterDestroy(scatter,ierr)
1534: !
1535: !  plot rval vector
1536: !
1537:       call VecView(rval,view0,ierr)
1538: !
1539: !  make xmgr plot of rval
1540: !
1541:       call FormXmgr(rval,1,ierr)
1542: !
1543: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1544: !
1545: !
1546: !  load eval from x
1547: !
1548:       call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1549:       call VecScatterBegin(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1550:       call VecScatterEnd(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1551:       call VecScatterDestroy(scatter,ierr)
1552: !
1553: !  plot eval vector
1554: !
1555:       call VecView(eval,view2,ierr)
1556: !
1557: !  make xmgr plot of eval
1558: !
1559:       call FormXmgr(eval,3,ierr)
1560: !
1561: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1562: !

1564: !
1565: !  load ruval from x
1566: !
1567:       call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1568:       call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1569:       call VecScatterEnd(scatter,x,ruval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1570:       call VecScatterDestroy(scatter,ierr)
1571: !
1572: !  create u = ru / r
1573: !
1574:       call VecPointwiseDivide(uval,ruval,rval,ierr)
1575: !
1576: !  plot uval vector
1577: !
1578:       call VecView(uval,view1,ierr)
1579: !
1580: !  make xmgr plot of uval
1581: !
1582:       call FormXmgr(uval,2,ierr)

1584: !
1585: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1586: !

1588:       call VecPointwiseMult(kval,uval,uval,ierr)
1589:       call VecScale(kval,0.5d+0,ierr)

1591:       call VecPointwiseDivide(steval,eval,rval,ierr)
1592:       call VecWAXPY(seval,-1.0d+0,kval,steval,ierr)

1594:       call VecCopy(seval,tval,ierr)
1595:       call VecScale(tval,csubvi,ierr)

1597: !
1598: !  plot tval vector
1599: !
1600:       call VecView(tval,view3,ierr)
1601: !
1602: !  make xmgr plot of tval
1603: !
1604:       call FormXmgr(tval,4,ierr)

1606: !
1607: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1608: !

1610:       call VecPointwiseMult(pval,rval,seval,ierr)
1611:       call VecScale(pval,gm1,ierr)
1612: !
1613: !  plot pval vector
1614: !
1615:       call VecView(pval,view4,ierr)
1616: !
1617: !  make xmgr plot of pval
1618: !
1619:       call FormXmgr(pval,5,ierr)
1620: !
1621: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1622: !





1628: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1629: !  Free work space.  All PETSc objects should be destroyed when they
1630: !  are no longer needed.
1631: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1633:       call VecDestroy(rval, ierr)
1634:       call VecDestroy(uval, ierr)
1635:       call VecDestroy(ruval,ierr)
1636:       call VecDestroy(eval, ierr)
1637:       call VecDestroy(seval, ierr)
1638:       call VecDestroy(pval, ierr)
1639:       call VecDestroy(kval, ierr)
1640:       call VecDestroy(tval, ierr)
1641:       call VecDestroy(steval, ierr)

1643:       call ISDestroy(rfrom, ierr)
1644:       call ISDestroy(rto,   ierr)

1646:       call ISDestroy(rufrom,ierr)
1647:       call ISDestroy(ruto,  ierr)

1649:       call ISDestroy(efrom, ierr)
1650:       call ISDestroy(eto,   ierr)


1653:       return
1654:       end
1655:       subroutine FormInitialGuess(x,ierr)
1656: ! ---------------------------------------------------------------------
1657: !
1658: !  FormInitialGuess - Forms initial approximation.
1659: !
1660: !  Input Parameter:
1661: !  x - vector
1662: !
1663: !  Output Parameters:
1664: !  x - vector
1665: !  ierr - error code
1666: !
1667: !  Notes:
1668: !  This routine serves as a wrapper for the lower-level routine
1669: !  "ApplicationInitialGuess", where the actual computations are
1670: !  done using the standard Fortran style of treating the local
1671: !  vector data as a multidimensional array over the local mesh.
1672: !  This routine merely accesses the local vector data via
1673: !  VecGetArray() and VecRestoreArray().
1674: !
1675:       implicit none

1677:  #include finclude/petsc.h

1679: !  Input/output variables:
1680:       Vec      x
1681:       integer  ierr

1683: !  Declarations for use with local arrays:
1684:       PetscScalar      lx_v(0:1)
1685:       PetscOffset lx_i

1687:       0

1689: !
1690: !  get a pointer to x
1691: !
1692:       call VecGetArray(x,lx_v,lx_i,ierr)
1693: !
1694: !  Compute initial guess
1695: !
1696:       call ApplicationInitialGuess(lx_v(lx_i),ierr)
1697: !
1698: !  Restore vector x
1699: !
1700:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1702:       return
1703:       end
1704:       subroutine FormXmgr(x,ivar,ierr)
1705: ! ---------------------------------------------------------------------
1706: !
1707: !  FormXmgr - Forms Xmgr output
1708: !
1709: !  Input Parameter:
1710: !  x - vector
1711: !
1712: !  Output Parameters:
1713: !  x - vector
1714: !  ierr - error code
1715: !
1716: !  Notes:
1717: !  This routine serves as a wrapper for the lower-level routine
1718: !  "ApplicationXmgr", where the actual computations are
1719: !  done using the standard Fortran style of treating the local
1720: !  vector data as a multidimensional array over the local mesh.
1721: !  This routine merely accesses the local vector data via
1722: !  VecGetArray() and VecRestoreArray().
1723: !
1724:       implicit none

1726:  #include finclude/petsc.h

1728: !  Input/output variables:
1729:       Vec      x
1730:       integer  ivar,ierr

1732: !  Declarations for use with local arrays:
1733:       PetscScalar      lx_v(0:1)
1734:       PetscOffset lx_i

1736:       0

1738: !
1739: !  get a pointer to x
1740: !
1741:       call VecGetArray(x,lx_v,lx_i,ierr)
1742: !
1743: !  make the graph
1744: !
1745:       call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1746: !
1747: !  Restore vector x
1748: !
1749:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1751:       return
1752:       end
1753:       subroutine PCRadApply(pc,ctx,x,y,ierr)
1754: ! -------------------------------------------------------------------
1755: !
1756: !   PCRadApply - This routine demonstrates the use of a
1757: !   user-provided preconditioner.
1758: !
1759: !   Input Parameters:
1760: !   dummy - optional user-defined context, not used here
1761: !   x - input vector
1762: !  In this example the shell preconditioner application context
1763: !  is a Fortran integer array:
1764: !      ctx(1) = shell preconditioner pressure matrix contex
1765: !      ctx(2) = semi implicit pressure matrix
1766: !      ctx(4) = xold  - old time values need for time advancement
1767: !      ctx(5) = mx    - number of control volumes
1768: !      ctx(6) = N     - total number of unknowns
1769: !
1770: !   Output Parameters:
1771: !   y - preconditioned vector
1772: !   ierr  - error code (nonzero if error has been detected)
1773: !
1774: !   Notes:
1775: !   This code implements the Jacobi preconditioner plus the
1776: !   SOR preconditioner
1777: !

1779:       implicit none

1781:  #include finclude/petsc.h

1783: #include "ex74fcomd.h"
1784: !
1785: !  Input
1786: !
1787:       PC               pc
1788:       PetscFortranAddr ctx(*)
1789:       Vec              x, y
1790:       integer          ierr
1791: !
1792: !  Local
1793: !
1794:       IS               defrom, deto
1795:       Vec              de, rese
1796:       VecScatter       scatter
1797:       PetscScalar           lde_v(0:1),lrese_v(0:1)
1798:       PetscOffset      lde_i,     lrese_i
1799: !
1800: !  Identity preconditioner
1801: !
1802:       call VecCopy(x,y,ierr)
1803: !
1804: !  if kappa0 not equal to zero then precondition the radiation diffusion
1805: !
1806:       if (kappa0 .ne. 0.0d+0) then
1807: 

1809: !
1810: !  Create needed vectors
1811: !
1812:          call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,de,ierr)
1813:          call VecSetFromOptions(de,ierr)
1814:          call VecDuplicate(de,rese,ierr)
1815: !
1816: !  create index sets for scatters
1817: !
1818:          call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1819:          call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1820: !
1821: !  load rese from x
1822: !
1823:          call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1824:          call VecScatterBegin(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1825:          call VecScatterEnd(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1826:          call VecScatterDestroy(scatter,ierr)
1827: !
1828: !  apply preconditioner
1829: !
1830:       call PCApply(ctx(1),rese,de,ierr)

1832:       if (debug) then
1833:         write(*,*) 'PCRadApply dh is'
1834:         call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1835:       endif
1836: !
1837: ! load de into y
1838: !
1839:       call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1840:       call VecScatterBegin(scatter,de,y,INSERT_VALUES, SCATTER_FORWARD,ierr)
1841:       call VecScatterEnd(scatter,de,y,INSERT_VALUES,SCATTER_FORWARD,ierr)
1842:       call VecScatterDestroy(scatter,ierr)

1844:       if (debug) then
1845:         write(*,*) 'PCRadApply y is'
1846:         call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1847:       endif

1849:       call VecDestroy(de,ierr)
1850:       call VecDestroy(rese,ierr)

1852:       call ISDestroy(defrom,ierr)
1853:       call ISDestroy(deto,ierr)

1855:       endif


1858:       return
1859:       end
1860:       subroutine PCRadSetUp(pc,ctx,ierr)
1861: !
1862: !   PCRadSetUp - This routine sets up a user-defined
1863: !   preconditioner context.
1864: !
1865: !   Input Parameters:
1866: !  In this example the shell preconditioner application context
1867: !  is a Fortran integer array:
1868: !      ctx(1) = shell preconditioner pressure matrix contex
1869: !      ctx(2) = semi implicit pressure matrix
1870: !      ctx(4) = xold  - old time values need for time advancement
1871: !      ctx(5) = mx    - number of control volumes
1872: !      ctx(6) = N     - total number of unknowns
1873: !
1874: !   Output Parameter:
1875: !   ierr  - error code (nonzero if error has been detected)
1876: !
1877: !   Notes:
1878: !   In this example, we define the shell preconditioner to be Jacobi
1879: !   method.  Thus, here we create a work vector for storing the reciprocal
1880: !   of the diagonal of the preconditioner matrix; this vector is then
1881: !   used within the routine PCRadApply().
1882: !

1884:       implicit none

1886:  #include finclude/petsc.h

1888: #include "ex74fcomd.h"
1889: !
1890: !  Input
1891: !
1892:       PC               pc
1893:       PetscFortranAddr ctx(*)
1894:       integer          ierr
1895: !
1896: !  Local
1897: !
1898:       Vec              eold
1899: 
1900:       PetscScalar      le_v(0:1)
1901:       PetscOffset le_i
1902: 
1903: !
1904: !  create vector
1905: !
1906:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,eold,ierr)
1907:       call VecSetFromOptions(eold,ierr)
1908: !
1909: ! set up the matrix based on xold
1910: !
1911:       call Setmat(ctx,ierr)
1912: !
1913: !  set up the preconditioner
1914: !
1915:       call PCDestroy(ctx(1),ierr)
1916:       call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
1917: !VAM  call PCSetType(ctx(1),PCJACOBI,ierr)
1918:       call PCSetType(ctx(1),PCLU,ierr)
1919: !      call PCSetVector(ctx(1),eold,ierr)
1920:       call PCSetOperators(ctx(1),ctx(2),ctx(2), ierr)
1921:       call PCSetUp(ctx(1),ierr)

1923:       call VecDestroy(eold,ierr)


1926:       return
1927:       end
1928:       subroutine Setmat(ctx,ierr)

1930:       implicit none

1932:  #include finclude/petsc.h

1934: !  Common blocks:
1935: #include "ex74fcomd.h"
1936: #include "ex74ftube.h"

1938: !  Input/output variables:
1939:       PetscFortranAddr ctx(*)
1940:       integer  ierr

1942: !  Local variables:
1943:       PetscScalar      lx_v(0:1)
1944:       PetscOffset lx_i

1946:       double precision xmult, himh, hiph, diag, upper, lower
1947:       double precision hi, hip1, him1
1948:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, ue,    up,    uw
1949:       double precision see, sep, sew, seef, sewf, tef, twf, ref, rwf, kef, kwf, xmulte, xmultw
1950: !
1951:       integer  im, nx, ny
1952: !
1953: !     get pointers to xold
1954: !
1955:       call VecGetArray(ctx(4),lx_v,lx_i,ierr)
1956: 

1958: !
1959: !############################
1960: !
1961: ! loop over all cells begin
1962: !
1963: !############################
1964: !
1965:       do im = 1,mx
1966: !
1967: !  set scalars
1968: !
1969:          call Setpbc(im,lx_v(lx_i), rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, ue,    up,    uw)
1970: !
1971: !  set diffusion coefficients
1972: !
1973:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
1974:         sep = (ergp/rhop) - (0.5d+0 * up * up)
1975:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

1977:         seef = 0.5d+0 * (see + sep)
1978:         sewf = 0.5d+0 * (sew + sep)

1980:         tef = seef / csubv
1981:         twf = sewf / csubv

1983:         ref = 0.5d+0 * (rhoe + rhop)
1984:         rwf = 0.5d+0 * (rhow + rhop)

1986:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
1987:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

1989:         if (wilson) then
1990:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
1991:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
1992:         endif
1993: !
1994: !  set coefficients
1995: !
1996:          xmult = dt / (dx * dx * csubv)

1998:          xmulte = xmult * kef
1999:          xmultw = xmult * kwf

2001:          upper = -(xmulte / rhoe)
2002:          lower = -(xmultw / rhow)

2004:          diag = 1.0d+0 + ( (xmulte + xmultw) / rhop )

2006: !
2007: !  load coefficients into the matrix
2008: !
2009:          call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)

2011:          if (im .eq. 1) then
2012:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper, INSERT_VALUES,ierr)
2013:          elseif (im .eq. mx) then
2014:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower,INSERT_VALUES,ierr)
2015:          else
2016:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper,INSERT_VALUES,ierr)
2017:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower, INSERT_VALUES,ierr)
2018:          endif


2021:       enddo
2022: !
2023: !############################
2024: !
2025: ! loop over all cells end
2026: !
2027: !############################
2028: !
2029: 
2030: !
2031: !  final load of matrix
2032: !
2033:       call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2034:       call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)

2036:       if (debug) then
2037:         call MatGetSize(ctx(2),nx,ny,ierr)
2038:         write(*,*) 'in setup nx = ',nx,' ny = ',ny
2039:         call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2040:       endif

2042:       call VecRestoreArray (ctx(4),lx_v,lx_i,ierr)



2046:       return
2047:       end
2048:       subroutine Setpbc(i,x, rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw)

2050:       implicit none

2052: !  Common blocks:
2053: #include "ex74fcomd.h"

2055: !  Input/output variables:
2056:       PetscScalar   x(mx*neq)
2057:       integer  i
2058:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2059:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2060:       double precision ergee,  erge,  ergp,  ergw,  ergww
2061:       double precision         vele,  velp,  velw

2063: !  Local variables:
2064:       integer  jr, jru, je

2066: !
2067: !  set pointers
2068: !
2069:       jr  = (neq*i) - 2
2070:       jru = (neq*i) - 1
2071:       je  = (neq*i)

2073:       if (debug) then
2074:         write(*,*)
2075:         write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2076:         write(*,*)
2077:       endif
2078: 
2079:       if (i .eq. 1) then

2081:         rhoee = x(jr+(2*neq))
2082:         rhoe  = x(jr+neq)
2083:         rhop  = x(jr)
2084:         rhow  = x(jr)
2085:         rhoww = x(jr)

2087:         rhouee = x(jru+(2*neq))
2088:         rhoue  = x(jru+neq)
2089:         rhoup  = x(jru)
2090:         rhouw  = x(jru)
2091:         rhouww = x(jru)

2093:         ergee = x(je+(2*neq))
2094:         erge  = x(je+neq)
2095:         ergp  = x(je)
2096:         ergw  = x(je)
2097:         ergww = x(je)

2099:         velw = 0.0d+0
2100:         velp = rhoup/rhop
2101:         vele = rhoue/rhoe

2103:       elseif (i .eq. 2) then

2105:         rhoee = x(jr+(2*neq))
2106:         rhoe  = x(jr+neq)
2107:         rhop  = x(jr)
2108:         rhow  = x(jr-neq)
2109:         rhoww = x(jr-neq)

2111:         rhouee = x(jru+(2*neq))
2112:         rhoue  = x(jru+neq)
2113:         rhoup  = x(jru)
2114:         rhouw  = x(jru-neq)
2115:         rhouww = x(jru-neq)

2117:         ergee = x(je+(2*neq))
2118:         erge  = x(je+neq)
2119:         ergp  = x(je)
2120:         ergw  = x(je-neq)
2121:         ergww = x(je-neq)

2123:         velw = rhouw/rhow
2124:         velp = rhoup/rhop
2125:         vele = rhoue/rhoe

2127:       elseif (i .eq. mx-1) then

2129:         rhoee = x(jr+neq)
2130:         rhoe  = x(jr+neq)
2131:         rhop  = x(jr)
2132:         rhow  = x(jr-neq)
2133:         rhoww = x(jr-(2*neq))

2135:         rhouee = x(jru+neq)
2136:         rhoue  = x(jru+neq)
2137:         rhoup  = x(jru)
2138:         rhouw  = x(jru-neq)
2139:         rhouww = x(jru-(2*neq))

2141:         ergee = x(je+neq)
2142:         erge  = x(je+neq)
2143:         ergp  = x(je)
2144:         ergw  = x(je-neq)
2145:         ergww = x(je-(2*neq))

2147:         velw = rhouw/rhow
2148:         velp = rhoup/rhop
2149:         vele = rhoue/rhoe

2151:       elseif (i .eq. mx) then

2153:         rhoee = x(jr)
2154:         rhoe  = x(jr)
2155:         rhop  = x(jr)
2156:         rhow  = x(jr-neq)
2157:         rhoww = x(jr-(2*neq))

2159:         rhouee = x(jru)
2160:         rhoue  = x(jru)
2161:         rhoup  = x(jru)
2162:         rhouw  = x(jru-neq)
2163:         rhouww = x(jru-(2*neq))

2165:         ergee = x(je)
2166:         erge  = x(je)
2167:         ergp  = x(je)
2168:         ergw  = x(je-neq)
2169:         ergww = x(je-(2*neq))

2171:         velw = rhouw/rhow
2172:         velp = rhoup/rhop
2173:         vele = 0.0d+0

2175:       else

2177:         rhoee = x(jr+(2*neq))
2178:         rhoe  = x(jr+neq)
2179:         rhop  = x(jr)
2180:         rhow  = x(jr-neq)
2181:         rhoww = x(jr-(2*neq))

2183:         rhouee = x(jru+(2*neq))
2184:         rhoue  = x(jru+neq)
2185:         rhoup  = x(jru)
2186:         rhouw  = x(jru-neq)
2187:         rhouww = x(jru-(2*neq))

2189:         ergee = x(je+(2*neq))
2190:         erge  = x(je+neq)
2191:         ergp  = x(je)
2192:         ergw  = x(je-neq)
2193:         ergww = x(je-(2*neq))

2195:         velw = rhouw/rhow
2196:         velp = rhoup/rhop
2197:         vele = rhoue/rhoe

2199:       endif

2201:       if (debug) then
2202:          write(*,*)
2203:          write(*,*) 'in Setpbc ',i,jr,jru,je
2204:          write(*,*) 'mx = ',mx
2205:          write(*,*)
2206:       endif


2209:       return
2210:       end
2211:       subroutine Setpbcn(rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, ue,    up,    uw,           jbc)

2213:       implicit none

2215: !  Common blocks:
2216: #include "ex74fcomd.h"

2218: !  Input/output variables:
2219:       integer  jbc
2220:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2221:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2222:       double precision ergee,  erge,  ergp,  ergw,  ergww
2223:       double precision         ue,    up,    uw

2225: !  Local variables:

2227:       if (jbc .eq. 1) then
2228:          rhoww  = rhop
2229:          rhow   = rhop
2230:          rhouww = rhoup
2231:          rhouw  = rhoup
2232:          ergww  = ergp
2233:          ergw   = ergp
2234:          uw     = 0.0d+0
2235:       elseif  (jbc .eq. 2) then
2236:          rhoww  = rhow
2237:          rhouww = rhouw
2238:          ergww  = ergw
2239:          uw     = rhouw / rhow
2240:       else
2241:          uw = rhouw / rhow
2242:       endif

2244:       if (jbc .eq. mx) then
2245:          rhoee  = rhop
2246:          rhoe   = rhop
2247:          rhouee = rhoup
2248:          rhoue  = rhoup
2249:          ergee  = ergp
2250:          erge   = ergp
2251:          ue     = 0.0d+0
2252:       elseif (jbc .eq. mx-1) then
2253:          rhoee  = rhoe
2254:          rhouee = rhoue
2255:          ergee  = erge
2256:          ue     = rhoue / rhoe
2257:       else
2258:          ue     = rhoue / rhoe
2259:       endif

2261:       up = rhoup / rhop

2263:       if (debug) then
2264:          write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2265:       endif


2268:       return
2269:       end
2270:       double precision function cont(rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww,jcont,xold)
2271: !
2272: !  This function computes the residual
2273: !  for the 1-D continuity equation
2274: !
2275: !
2276:       implicit none

2278:       include 'ex74fcomd.h'
2279:       include 'ex74ftube.h'
2280: !
2281: !     input variables
2282: !
2283:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2284:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2285:       double precision ergee,  erge,  ergp,  ergw,  ergww
2286:       double precision xold(mx*neq)
2287: !
2288:       integer jcont
2289: !
2290: !     local variables
2291: !
2292:       double precision theta1
2293:       integer jr
2294: !
2295: !  new
2296: !
2297:       double precision velfw, velfe
2298:       double precision vele,velp,velw
2299:       double precision fluxe, fluxw
2300:       double precision urhoe, urhow
2301:       double precision source
2302: !
2303: ! old
2304: !
2305:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2306:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2307:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2308:       double precision teoee, teoe, teop, teow, teoww, uoe, uoee, uop, uow, uoww
2309:       double precision velfow, velfoe
2310:       double precision veloe,velop,velow
2311:       double precision fluxoe, fluxow
2312:       double precision urhooe, urhoow
2313:       double precision sourceo
2314: !
2315: ! functions
2316: !
2317:       double precision godunov2
2318:       double precision upwind, fluxlim
2319: !
2320: !
2321: ! ******************************************************************
2322: !
2323: !
2324: !
2325:       if (debug) then
2326:         write(*,*)
2327:         write(*,*) 'in cont',jcont,' ihod = ',ihod
2328:         write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2329:         write(*,*) 'rhop = ',rhop
2330:         write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2331:         write(*,*)
2332:       endif

2334:       jr = (neq*jcont) - 2

2336: !########################
2337: !
2338: !      NEW
2339: !
2340: !########################

2342:       call Setpbcn(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw,         jcont)

2344:       velfe = 0.5d+0 * (vele + velp)
2345:       velfw = 0.5d+0 * (velw + velp)

2347:       if (ihod .eq. 1) then

2349:         urhoe = upwind(rhop,rhoe,velfe)
2350:         urhow = upwind(rhow,rhop,velfw)

2352:       elseif (ihod .eq. 2) then

2354:         urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2355:         urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)

2357:       endif

2359:       if (ihod .eq. 3) then
2360:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,1)
2361:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,1)
2362:       else
2363:         fluxe = (dt/dx) * urhoe
2364:         fluxw = (dt/dx) * urhow
2365:       endif


2368:       source = 0.0d+0

2370: !########################
2371: !
2372: !      OLD
2373: !
2374: !########################

2376:       call Setpbc(jcont,xold,rhooee,  rhooe,  rhoop,  rhoow,  rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow)

2378:       call Setpbcn(rhooee,  rhooe,  rhoop,  rhoow,  rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow,         jcont)

2380:       velfoe = 0.5d+0 * (veloe + velop)
2381:       velfow = 0.5d+0 * (velow + velop)


2384:       if (ihod .eq. 1) then

2386:         urhooe = upwind(rhoop,rhooe,velfoe)
2387:         urhoow = upwind(rhoow,rhoop,velfow)

2389:       elseif (ihod .eq. 2) then

2391:         urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2392:         urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)

2394:       endif

2396:       if (ihod .eq. 3) then
2397:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,1)
2398:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,1)
2399:       else
2400:         fluxoe = (dt/dx) * urhooe
2401:         fluxow = (dt/dx) * urhoow
2402:       endif

2404:       sourceo = 0.0d+0


2407: !########################
2408: !
2409: !      FUNCTION
2410: !
2411: !########################

2413:       theta1 = 1.0d+0 - theta
2414:       cont =  (rhop - xold(jr))  + (  theta  * ( (fluxe  - fluxw ) - source  )  )   + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
2415: !VAM
2416:       if (probnum .eq. 3) then
2417:         cont = 0.0d+0
2418:       endif
2419: !VAM


2422:       if (debug) then
2423:        write(*,*)
2424:        write(*,*) 'cont(',jcont,') = ',cont
2425:        write(*,*) 'theta = ',theta,'rhop = ',rhop
2426:        write(*,*) 'source = ',source,' sourceo = ',sourceo
2427:        write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2428:        write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2429:        write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2430:        write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2431:        write(*,*)
2432:       endif

2434:       return
2435:       end
2436:       double precision function  eexact(x,t)

2438:       implicit none

2440:       double precision x,t
2441:       double precision xot, head, tail, contact, ufan
2442:       double precision xpow, grat, urat
2443:       double precision uexact


2446:       logical debug

2448:       include 'ex74ftube.h'

2450:       debug = .false.


2453:       if (t .le. 0.0d+0) then
2454:         if (x .gt. 0.0d+0) then
2455:           eexact = e1
2456:         else
2457:           eexact = e4
2458:         endif
2459:       else

2461:        xot = x/t
2462:        head = -a4
2463:        tail = v3 - a3
2464:        contact = v2

2466:        if (xot .lt. head) then
2467:           eexact = e4
2468:        elseif (xot .gt. sspd) then
2469:           eexact = e1
2470:        elseif (xot .gt. contact) then
2471:           eexact = e2
2472:        elseif (xot .gt. tail) then
2473:           eexact = e3
2474:        else
2475:           ufan = uexact(x,t)
2476:           grat = (gamma - 1.0d+0) / 2.0d+0
2477:           xpow = 2.0d+0
2478:           urat = ufan / a4
2479:           eexact = e4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
2480:        endif

2482:       endif


2485:       if (debug) then
2486:         write(*,*)
2487:         write(*,*) 'eexact(',x,',',t,') = ',eexact
2488:         write(*,*)
2489:       endif

2491:       return
2492:       end
2493:       subroutine eigen(ht,uht)
2494: !23456789012345678901234567890123456789012345678901234567890123456789012
2495: !
2496: !          subroutine eigen
2497: !
2498: !  This subroutine computes the eigen values and eigen vectors
2499: !
2500: !23456789012345678901234567890123456789012345678901234567890123456789012


2503: !#######################################################################

2505:       implicit none

2507:       include 'ex74fcomd.h'

2509:       double precision ht, uht

2511:       double precision ut, at, lam1, lam2


2514: !#######################################################################

2516:       ut = uht / ht
2517:       at = sqrt( ht)

2519:       lam1 = ut - at
2520:       lam2 = ut + at

2522:       eigval(1) = lam1
2523:       eigval(2) = lam2

2525:       eigvec(1,1) = 1.0d+0
2526:       eigvec(2,1) = lam1
2527:       eigvec(1,2) = 1.0d+0
2528:       eigvec(2,2) = lam2

2530:       rinv(1,1) =  lam2 / (2.0d+0 * at)
2531:       rinv(2,1) = -lam1 / (2.0d+0 * at)
2532:       rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2533:       rinv(2,2) =  1.0d+0 / (2.0d+0 * at)


2536:       return
2537:       end
2538:       subroutine eigene(r,ru,e,l1,l2,l3)
2539: !23456789012345678901234567890123456789012345678901234567890123456789012
2540: !
2541: !          subroutine eigene
2542: !
2543: !  This subroutine computes the eigen values for the entropy fix
2544: !
2545: !23456789012345678901234567890123456789012345678901234567890123456789012


2548: !#######################################################################

2550:       implicit none

2552:       include 'ex74ftube.h'

2554:       double precision r,ru,e,l1,l2,l3

2556:       double precision p,u,a

2558:       double precision eos
2559:       integer ierr

2561:       logical debug


2564: !#######################################################################

2566:       debug = .false.

2568:       if (debug) then
2569:          write(*,*)
2570:          write(*,*) 'gamma = ',gamma
2571:          write(*,*) 'r,ru,e = ',r,ru,e
2572:          write(*,*)
2573:       endif

2575:       p = eos(r,ru,e)
2576:       u = ru/r
2577:       if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2578:          write(*,*)
2579:          write(*,*) 'gamma = ',gamma
2580:          write(*,*) 'r = ',r
2581:          write(*,*) 'p = ',p
2582:          write(*,*)
2583:          call PetscFinalize(ierr)
2584:          stop
2585:       endif
2586:       a = sqrt((gamma * p)/r)

2588:       if (debug) then
2589:          write(*,*)
2590:          write(*,*) 'p,u,a = ',p,u,a
2591:          write(*,*)
2592:       endif

2594:       l1 = u - a
2595:       l2 = u
2596:       l3 = u + a

2598:       return
2599:       end
2600:       double precision function energy(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,jerg,xold)
2601: !
2602: !  This function computes the residual
2603: !  for the 1-D energy equation
2604: !
2605: !
2606:       implicit none

2608:       include 'ex74fcomd.h'
2609:       include 'ex74ftube.h'
2610: !
2611: !     input variables
2612: !
2613:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2614:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2615:       double precision ergee,  erge,  ergp,  ergw,  ergww
2616:       double precision xold(mx*neq)
2617: !
2618:       integer jerg
2619: !
2620: !     local variables
2621: !
2622:       double precision theta1
2623:       integer je
2624: !
2625: !  new
2626: !
2627:       double precision velfw, velfe
2628:       double precision vele,velp,velw
2629:       double precision fluxe, fluxw
2630:       double precision uepe, uepw
2631:       double precision ue, up, uw
2632:       double precision see, sep, sew
2633:       double precision seef, sewf
2634:       double precision upe, upw
2635:       double precision presse, pressw
2636:       double precision source
2637:       double precision te, tp, tw
2638:       double precision tef, twf, ref, rwf
2639:       double precision kef, kwf
2640:       double precision hflxe, hflxw
2641: !
2642: ! old
2643: !
2644:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2645:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2646:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2647:       double precision velfow, velfoe
2648:       double precision veloe,velop,velow
2649:       double precision fluxoe, fluxow
2650:       double precision uepoe, uepow
2651:       double precision uoe, uop, uow
2652:       double precision seoe, seop, seow
2653:       double precision seoef, seowf
2654:       double precision upoe, upow
2655:       double precision pressoe, pressow
2656:       double precision sourceo
2657:       double precision toe, top, tow
2658:       double precision toef, towf, roef, rowf
2659:       double precision koef, kowf
2660:       double precision hflxoe, hflxow
2661: !
2662: ! functions
2663: !
2664:       double precision godunov2, eos
2665:       double precision upwind, fluxlim

2667: !
2668: !
2669: ! ******************************************************************
2670: !
2671: !
2672: !
2673:       je = (neq*jerg)

2675: !########################
2676: !
2677: !      NEW
2678: !
2679: !########################

2681:       call Setpbcn(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,vele,  velp,  velw,         jerg)

2683:       pressw  = eos(rhow, rhouw, ergw)
2684:       presse  = eos(rhoe, rhoue, erge)

2686:       uw = rhouw / rhow
2687:       up = rhoup / rhop
2688:       ue = rhoue / rhoe

2690:       upw = uw * pressw
2691:       upe = ue * presse

2693:       velfe = 0.5d+0 * (vele + velp)
2694:       velfw = 0.5d+0 * (velw + velp)

2696:       if (ihod .eq. 1) then

2698:         uepe = upwind(ergp,erge,velfe)
2699:         uepw = upwind(ergw,ergp,velfw)

2701:       elseif (ihod .eq. 2) then

2703:         uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2704:         uepw = fluxlim(ergww,ergw,ergp,erge,velfw)

2706:       endif

2708:       if (ihod .eq. 3) then
2709:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,3)
2710:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,3)
2711:       else
2712:         fluxe = (dt/dx) * ( uepe  + (0.5d+0*upe) )
2713:         fluxw = (dt/dx) * ( uepw  + (0.5d+0*upw) )
2714:       endif
2715: !
2716: !  radiation
2717: !
2718:       if (kappa0 .eq. 0.0d+0) then
2719:         source = 0.0d+0
2720:       else

2722:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
2723:         sep = (ergp/rhop) - (0.5d+0 * up * up)
2724:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

2726:         seef = 0.5d+0 * (see + sep)
2727:         sewf = 0.5d+0 * (sew + sep)

2729:         te  = see / csubv
2730:         tp  = sep / csubv
2731:         tw  = sew / csubv

2733:         tef = seef / csubv
2734:         twf = sewf / csubv

2736:         ref = 0.5d+0 * (rhoe + rhop)
2737:         rwf = 0.5d+0 * (rhow + rhop)

2739:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2740:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

2742:         if (wilson) then
2743:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2744:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2745:         endif

2747:         if ( debug .and. (kef .gt. 1.0d+10) ) then
2748:           write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2749:         endif
2750:         if ( debug .and. (kwf .gt. 1.0d+10) ) then
2751:           write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2752:         endif

2754:         hflxe = kef * (te - tp) / dx
2755:         hflxw = kwf * (tp - tw) / dx

2757:         source = (dt/dx) * (hflxe - hflxw)

2759:       endif

2761: !########################
2762: !
2763: !      OLD
2764: !
2765: !########################

2767:       call Setpbc(jerg,xold, rhooee,  rhooe,  rhoop,  rhoow,  rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow)

2769:       call Setpbcn(rhooee,  rhooe,  rhoop,  rhoow,  rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow,         jerg)

2771:       pressow  = eos(rhoow, rhouow, ergow)
2772:       pressoe  = eos(rhooe, rhouoe, ergoe)

2774:       uow = rhouow / rhoow
2775:       uop = rhouop / rhoop
2776:       uoe = rhouoe / rhooe

2778:       upow = uow * pressow
2779:       upoe = uoe * pressoe

2781:       velfoe = 0.5d+0 * (veloe + velop)
2782:       velfow = 0.5d+0 * (velow + velop)


2785:       if (ihod .eq. 1) then

2787:         uepoe = upwind(ergop,ergoe,velfoe)
2788:         uepow = upwind(ergow,ergop,velfow)

2790:       elseif (ihod .eq. 2) then

2792:         uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
2793:         uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)

2795:       endif

2797:       if (ihod .eq. 3) then
2798:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,3)
2799:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,3)
2800:       else
2801:         fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
2802:         fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
2803:       endif

2805: !
2806: !  old radiation
2807: !
2808:       if (kappa0 .eq. 0.0d+0) then
2809:         sourceo = 0.0d+0
2810:       else

2812:         seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
2813:         seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
2814:         seow = (ergow/rhoow) - (0.5d+0 * uow * uow)

2816:         seoef = 0.5d+0 * (seoe + seop)
2817:         seowf = 0.5d+0 * (seow + seop)

2819:         toe  = seoe / csubv
2820:         top  = seop / csubv
2821:         tow  = seow / csubv

2823:         toef = seoef / csubv
2824:         towf = seowf / csubv

2826:         roef = 0.5d+0 * (rhooe + rhoop)
2827:         rowf = 0.5d+0 * (rhoow + rhoop)

2829:         koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
2830:         kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)

2832:         if (wilson) then
2833:            koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
2834:            kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
2835:         endif

2837:         if ( debug .and. (koef .gt. 1.0d+10) ) then
2838:           write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
2839:         endif
2840:         if ( debug .and. (kowf .gt. 1.0d+10) ) then
2841:           write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
2842:         endif

2844:         hflxoe = koef * (toe - top) / dx
2845:         hflxow = kowf * (top - tow) / dx

2847:         sourceo = (dt/dx) * (hflxoe - hflxow)

2849:       endif


2852: !########################
2853: !
2854: !      FUNCTION
2855: !
2856: !########################

2858: !VAM
2859:       if (probnum .eq. 3) then
2860:         fluxe  = 0.0d+0
2861:         fluxw  = 0.0d+0
2862:         fluxoe = 0.0d+0
2863:         fluxow = 0.0d+0
2864:       endif
2865: !VAM

2867:       theta1 = 1.0d+0 - theta
2868:       energy =  (ergp - xold(je))  + (  theta  * ( (fluxe  - fluxw ) - source  )  )  + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )

2870:       if (debug) then
2871:         write(*,*)
2872:         write(*,*) 'energy(',jerg,') = ',energy
2873:         write(*,*)
2874:         write(*,*) fluxe,fluxw
2875:         write(*,*) fluxoe,fluxow
2876:         write(*,*) source,sourceo
2877:         write(*,*)
2878:       endif

2880:       return
2881:       end
2882:       double precision function eos(r,ru,e)

2884:       implicit none

2886:       double precision r,ru,e

2888:       double precision se, u

2890:       integer ierr

2892:       logical debug

2894:       include "ex74ftube.h"

2896:       debug = .false.

2898:       if (debug) then
2899:         write(*,*)
2900:         write(*,*) 'in eos r,ru,e'
2901:         write(*,*) r,ru,e
2902:         write(*,*)
2903:       endif

2905:       u = ru/r

2907:       se = (e/r) - (0.5d+0 * u * u)
2908:       eos = (gamma - 1.0d+0) * r * se

2910:       if (eos .lt. 0.0d+0) then
2911:          write(*,*)
2912:          write(*,*) 'eos = ',eos
2913:          write(*,*) 'gamma = ',gamma
2914:          write(*,*) 'r = ',r
2915:          write(*,*) 'se = ',se
2916:          write(*,*) 'e = ',e
2917:          write(*,*) 'u = ',u
2918:          write(*,*) 'ru = ',ru
2919:          call PetscFinalize(ierr)
2920:          write(*,*)
2921:          stop
2922:       endif

2924:       if (debug) then
2925:         write(*,*)
2926:         write(*,*) 'in eos u,se,eos'
2927:         write(*,*) u,se,eos
2928:         write(*,*)
2929:       endif


2932:       return
2933:       end
2934:       subroutine eval2

2936:       implicit none

2938:       double precision prat, grat, xnum, xdenom


2941:       logical debug

2943:       include 'ex74ftube.h'

2945:       debug = .false.

2947:       prat = p2/p1
2948:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

2950:       xnum = grat + prat
2951:       xdenom = 1.0d+0 + (prat * grat)
2952: 
2953:       e2 = e1 * prat * (xnum/xdenom)
2954: 


2957:       if (debug) then
2958:         write(*,*)
2959:         write(*,*) 'e1  = ',e1
2960:         write(*,*) 'e2  = ',e2
2961:       endif

2963:       return
2964:       end
2965:       subroutine exact0

2967:       implicit none

2969:       double precision tol, xn
2970:       double precision shockp, fprime

2972:       integer maxnewt, niter

2974:       logical found, debug

2976:       include 'ex74ftube.h'

2978:       debug = .false.

2980:       tol = 1.0d-10

2982:       maxnewt = 40
2983: 
2984:       a1 = sqrt(gamma*p1/r1)
2985:       a4 = sqrt(gamma*p4/r4)



2989:       found = .false.
2990:       niter = 0

2992:       xn =  0.5d+0 * (p1 + p4)

2994:    10 if ( (.not. found) .and. (niter .le. maxnewt) ) then

2996:         niter = niter + 1

2998:         xn = xn - (shockp(xn) / fprime(xn))

3000:         if (debug) then
3001:           write(*,*) niter,shockp(xn),xn
3002:         endif

3004:         if ( abs(shockp(xn)) .lt. tol ) then
3005:            found = .true.
3006:         endif

3008:         goto 10

3010:       endif

3012:       if (.not. found) then

3014:          write(*,*) 'newton failed'
3015:          write(*,*) xn,shockp(xn)
3016:          stop

3018:       endif

3020:       p2 = xn


3023:       if (debug) then
3024:         write(*,*)
3025:         write(*,*) 'p1  = ',p1
3026:         write(*,*) 'p2  = ',p2
3027:         write(*,*) 'p4  = ',p4
3028:         write(*,*)
3029:       endif

3031:       return
3032:       end
3033:       double precision function flux(r,ru,e,eqn)
3034: !23456789012345678901234567890123456789012345678901234567890123456789012
3035: !
3036: !          function flux
3037: !
3038: !  This function computes the flux at a face
3039: !
3040: !23456789012345678901234567890123456789012345678901234567890123456789012


3043: !#######################################################################

3045:       implicit none

3047:       include 'ex74fcomd.h'
3048:       include 'ex74ftube.h'

3050:       double precision r, ru, e

3052:       integer eqn

3054:       double precision p,u

3056:       double precision eos


3059: !#######################################################################

3061:       p = eos(r,ru,e)
3062:       u = ru/r

3064:       if (eqn .eq. 1) then
3065:          flux = ru
3066:       elseif (eqn .eq. 2) then
3067:          flux = (u * ru) + p
3068:       else
3069:          flux = u * (e + p)
3070:       endif

3072:       return
3073:       end
3074:       double precision function fluxlim(fww,fw,fe,fee,vp)
3075: !23456789012345678901234567890123456789012345678901234567890123456789012
3076: !
3077: !          function fluxlim
3078: !
3079: !  this function computes the flux limited quick face value
3080: !
3081: !23456789012345678901234567890123456789012345678901234567890123456789012


3084: !#######################################################################

3086:       implicit none

3088:       double precision fww, fw, fe, fee, vp

3090:       double precision fd, fc, fu

3092:       double precision f1, f2, f3, f4, fhod, beta, flc

3094:       double precision med, quick
3095: 
3096:       logical limit

3098: !#######################################################################

3100:       limit = .true.

3102:       if (vp .gt. 0.0d+0) then
3103:         fd = fe
3104:         fc = fw
3105:         fu = fww
3106:       else
3107:         fd = fw
3108:         fc = fe
3109:         fu = fee
3110:       endif

3112:       fhod = quick(fd,fc,fu)

3114:       if (limit) then

3116:         beta = 0.25d+0
3117:         flc = 4.0d+0

3119:         f1 = fc
3120:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3121:         f3 = fu + ( flc * (fc - fu) )
3122:         f4 = med(f1,f2,f3)
3123:         fluxlim = vp * med(f1,f4,fhod)

3125:       else

3127:         fluxlim = vp * fhod

3129:       endif

3131:       return
3132:       end
3133:       double precision function fluxlim2(fww,fw,fe,fee,vp)
3134: !23456789012345678901234567890123456789012345678901234567890123456789012
3135: !
3136: !          function fluxlim2
3137: !
3138: !  this function computes the flux limited quick face value
3139: !
3140: !23456789012345678901234567890123456789012345678901234567890123456789012


3143: !#######################################################################

3145:       implicit none

3147:       double precision fww, fw, fe, fee, vp

3149:       double precision fd, fc, fu

3151:       double precision f1, f2, f3, f4, fhod, beta, flc

3153:       double precision med, quick
3154: 
3155:       logical limit, debug

3157: !#######################################################################

3159:       debug = .false.

3161:       if (debug) then
3162:         write(*,*)
3163:         write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3164:         write(*,*) fee,fe,fw,fww
3165:         write(*,*)
3166:       endif

3168:       limit = .true.

3170:       if (vp .gt. 0.0d+0) then
3171:         fd = fe
3172:         fc = fw
3173:         fu = fww
3174:       else
3175:         fd = fw
3176:         fc = fe
3177:         fu = fee
3178:       endif

3180:       fhod = quick(fd,fc,fu)

3182:       if (limit) then

3184:         beta = 0.25d+0
3185:         flc = 4.0d+0

3187:         f1 = fc
3188:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3189:         f3 = fu + ( flc * (fc - fu) )
3190:         f4 = med(f1,f2,f3)
3191:         fluxlim2 =  med(f1,f4,fhod)

3193:       else

3195:         fluxlim2 = fhod

3197:       endif

3199:       return
3200:       end
3201:       double precision function fprime(x)

3203:       implicit none

3205:       double precision  x, eps
3206:       double precision  shockp

3208:       eps = 1.0d-8

3210:       fprime = ( shockp(x+eps) - shockp(x) ) / eps

3212:       return
3213:       end
3214:       double precision function godunov2(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3215: !23456789012345678901234567890123456789012345678901234567890123456789012
3216: !
3217: !          function godunov2
3218: !
3219: !  this function computes the roe/godunov2 face value
3220: !
3221: !23456789012345678901234567890123456789012345678901234567890123456789012


3224: !#######################################################################

3226:       implicit none

3228:       include 'ex74fcomd.h'
3229:       include 'ex74ftube.h'

3231:       double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err

3233:       integer eqn

3235:       double precision rrg, rlg, rurg, rulg, erg, elg

3237:       double precision hlle


3240: !#######################################################################

3242:       if (gorder .eq. 1) then
3243:         rrg  = rr
3244:         rlg  = rl
3245:         rurg = rur
3246:         rulg = rul
3247:         erg  = er
3248:         elg  = el
3249:       else
3250:         call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,rrg, rlg,rurg, rulg, erg, elg)
3251:       endif

3253: !VAM  if (ientro .eq. 0) then
3254: !VAM     godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3255: !VAM  elseif(ientro .eq. 1) then
3256: !VAM     godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3257: !VAM  else
3258:          godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3259: !VAM  endif


3262:       return
3263:       end
3264:       double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3265: !23456789012345678901234567890123456789012345678901234567890123456789012
3266: !
3267: !          function hlle
3268: !
3269: !  this function computes the roe/hlle face value
3270: !
3271: !23456789012345678901234567890123456789012345678901234567890123456789012


3274: !#######################################################################

3276:       implicit none

3278:       include 'ex74fcomd.h'
3279:       include 'ex74ftube.h'

3281:       double precision rrg,rlg,rurg,rulg,erg,elg
3282:       integer eqn

3284:       double precision laml1, laml2, laml3
3285:       double precision lamr1, lamr2, lamr3
3286:       double precision sl, sr


3289:       double precision flux

3291:       integer i, j, ispeed
3292: 

3294: !#######################################################################

3296:       ispeed = 1

3298:       do i = 1,neq
3299:         fr(i) = flux(rrg,rurg,erg,i)
3300:         fl(i) = flux(rlg,rulg,elg,i)
3301:       enddo

3303:       deltau(1) = rrg  - rlg
3304:       deltau(2) = rurg - rulg
3305:       deltau(3) = erg  - elg

3307: !VAM  call roestat(uhl,uhr,hl,hr,ht,uht)

3309: !VAM  call eigene(ht,uht,lamt1, lamt2)
3310:       call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3311:       call eigene(rlg,rulg,elg,laml1,laml2,laml3)

3313: !VAM  if (ispeed .eq. 1) then
3314: !VAM    sl = min(laml1,lamt1)
3315: !VAM    sr = max(lamt2,lamr2)
3316: !VAM  else
3317:         sl = min(laml1,lamr1)
3318:         sr = max(laml3,lamr3)
3319: !VAM  endif


3322:       do i = 1,neq
3323:         froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )/(sr-sl)
3324:       enddo

3326:       hlle = froe(eqn)


3329:       return
3330:       end
3331:       double precision function med(x1,x2,x3)
3332: !23456789012345678901234567890123456789012345678901234567890123456789012
3333: !
3334: !          function med
3335: !
3336: !  this function computes the median of three numbers
3337: !
3338: !23456789012345678901234567890123456789012345678901234567890123456789012


3341: !#######################################################################

3343:       implicit none

3345:       double precision x1, x2, x3
3346:       double precision xhi, xlo

3348: !#######################################################################

3350:       xhi = max(x1,x2,x3)
3351:       xlo = min(x1,x2,x3)

3353:       med = x1 + x2 + x3 - xhi - xlo

3355:       return
3356:       end
3357:       double precision function mom(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,jmom,xold)
3358: !
3359: !  This function computes the residual
3360: !  for the 1-D momentum equation
3361: !
3362: !
3363:       implicit none

3365:       include 'ex74fcomd.h'
3366:       include 'ex74ftube.h'
3367: !
3368: !     input variables
3369: !
3370:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
3371:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
3372:       double precision ergee,  erge,  ergp,  ergw,  ergww
3373:       double precision xold(mx*neq)
3374: !
3375:       integer jmom
3376: !
3377: !     local variables
3378: !
3379:       double precision theta1
3380:       integer jru
3381: !
3382: !  new
3383: !
3384:       double precision velfw, velfe
3385:       double precision vele,velp,velw
3386:       double precision fluxe, fluxw
3387:       double precision uurhoe, uurhow
3388:       double precision pressee, presse, pressp,pressw, pressww
3389:       double precision rupee, rupe, rupp, rupw, rupww
3390:       double precision uee, ue, up, uw, uww
3391:       double precision source
3392: !
3393: ! old
3394: !
3395:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
3396:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3397:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
3398:       double precision velfow, velfoe
3399:       double precision veloe,velop,velow
3400:       double precision fluxoe, fluxow
3401:       double precision uurhooe, uurhoow
3402:       double precision pressoee, pressoe, pressop, pressow, pressoww
3403:       double precision rupoee, rupoe, rupop, rupow, rupoww
3404:       double precision uoee, uoe, uop, uow, uoww
3405:       double precision sourceo

3407:       double precision eps
3408: !
3409: ! functions
3410: !
3411:       double precision godunov2, eos
3412:       double precision upwind, fluxlim
3413: !
3414: !
3415: ! ******************************************************************
3416: !
3417: !
3418:       eps = 1.0d-32
3419: !
3420:       jru = (neq*jmom) - 1

3422: !########################
3423: !
3424: !      NEW
3425: !
3426: !########################

3428:       call Setpbcn(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,vele,  velp,  velw, jmom)

3430:       presse  = eos(rhoe, rhoue, erge )
3431:       pressw  = eos(rhow, rhouw, ergw )

3433:       velfe = 0.5d+0 * (vele + velp)
3434:       velfw = 0.5d+0 * (velw + velp)

3436:       if (ihod .eq. 1) then

3438:         uurhoe = upwind(rhoup,rhoue,velfe)
3439:         uurhow = upwind(rhouw,rhoup,velfw)

3441:       elseif (ihod .eq. 2) then

3443:         uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3444:         uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)

3446:       endif

3448:       if (ihod .eq. 3) then
3449:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,2)
3450:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,2)
3451:       else
3452:         fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3453:         fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3454:       endif


3457:       source = 0.0d+0

3459: !########################
3460: !
3461: !      OLD
3462: !
3463: !########################

3465:       call Setpbc(jmom,xold,rhooee,  rhooe,  rhoop,  rhoow,  rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow)

3467:       call Setpbcn(rhooee,  rhooe,  rhoop,  rhoow,  rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee,  ergoe,  ergop,  ergow,  ergoww, veloe,  velop,  velow, jmom)

3469:       pressoe  = eos(rhooe, rhouoe, ergoe)
3470:       pressow  = eos(rhoow, rhouow, ergow)

3472:       velfoe = 0.5d+0 * (veloe + velop)
3473:       velfow = 0.5d+0 * (velow + velop)

3475:       if (ihod .eq. 1) then

3477:         uurhooe = upwind(rhouop,rhouoe,velfoe)
3478:         uurhoow = upwind(rhouow,rhouop,velfow)

3480:       elseif (ihod .eq. 2) then

3482:         uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3483:         uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)

3485:       endif

3487:       if (ihod .eq. 3) then
3488:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee,ergow, ergop, ergoe, ergoee,2)
3489:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,2)
3490:       else
3491:         fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3492:         fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3493:       endif

3495:       sourceo = 0.0d+0


3498: !########################
3499: !
3500: !      FUNCTION
3501: !
3502: !########################

3504:       theta1 = 1.0d+0 - theta
3505:       mom =  (rhoup - xold(jru))  + (  theta  * ( (fluxe  - fluxw ) - source  )  )  + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
3506: !VAM
3507:       if (probnum .eq. 3) then
3508:         mom = 0.0d+0
3509:       endif
3510: !VAM
3511:       if (debug) then
3512:         write(*,*)
3513:         write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3514:         write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3515:         write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3516:         write(*,*) 'presse = ',presse,'pressw = ',pressw
3517:         write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3518:         write(*,*)
3519:       endif

3521:       return
3522:       end
3523:       double precision function quick(fd, fc, fu)
3524: !23456789012345678901234567890123456789012345678901234567890123456789012
3525: !
3526: !          function quick
3527: !
3528: !  this function computes the quick face value
3529: !
3530: !23456789012345678901234567890123456789012345678901234567890123456789012


3533: !#######################################################################

3535:       implicit none

3537:       double precision fd, fc, fu

3539: !#######################################################################

3541:       quick = ( (3.0d+0 * fd) + (6.0d+0 * fc) - fu ) / 8.0d+0

3543:       return
3544:       end
3545:       double precision function  rexact(x,t)

3547:       implicit none

3549:       double precision x,t
3550:       double precision xot, head, tail, contact, ufan
3551:       double precision xpow, grat, urat
3552:       double precision uexact


3555:       logical debug

3557:       include 'ex74ftube.h'

3559:       debug = .false.


3562:       if (t .le. 0.0d+0) then
3563:         if (x .gt. 0.0d+0) then
3564:           rexact = r1
3565:         else
3566:           rexact = r4
3567:         endif
3568:       else

3570:        xot = x/t
3571:        head = -a4
3572:        tail = v3 - a3
3573:        contact = v2

3575:        if (xot .lt. head) then
3576:           rexact = r4
3577:        elseif (xot .gt. sspd) then
3578:           rexact = r1
3579:        elseif (xot .gt. contact) then
3580:           rexact = r2
3581:        elseif (xot .gt. tail) then
3582:           rexact = r3
3583:        else
3584:           ufan = uexact(x,t)
3585:           grat = (gamma - 1.0d+0) / 2.0d+0
3586:           xpow = 1.0d+0 / grat
3587:           urat = ufan / a4
3588:           rexact = r4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
3589:        endif

3591:       endif


3594:       if (debug) then
3595:         write(*,*)
3596:         write(*,*) 'rexact(',x,',',t,') = ',rexact
3597:         write(*,*)
3598:       endif

3600:       return
3601:       end
3602:       subroutine roestat(uhl, uhr, hl,hr,ht,uht)
3603: !23456789012345678901234567890123456789012345678901234567890123456789012
3604: !
3605: !          subroutine roestat
3606: !
3607: !  This subroutine computes the roe state at a face
3608: !
3609: !23456789012345678901234567890123456789012345678901234567890123456789012


3612: !#######################################################################

3614:       implicit none

3616:       include 'ex74fcomd.h'

3618:       double precision uhl, uhr, hl, hr, ht, uht

3620:       double precision ul, ur, shl, shr, xnum, xdenom
3621: 


3624: !#######################################################################

3626:       ul = uhl / hl
3627:       ur = uhr / hr

3629:       shl = sqrt(hl)
3630:       shr = sqrt(hr)

3632:       xnum = (shl * ul) + (shr * ur)
3633:       xdenom = shl + shr

3635:       ht  = 0.5d+0 * (hl + hr)
3636:       uht = ht * ( xnum / xdenom )

3638:       return
3639:       end
3640:       subroutine rval2

3642:       implicit none

3644:       double precision prat, grat, xnum, xdenom


3647:       logical debug

3649:       include 'ex74ftube.h'

3651:       debug = .false.

3653:       prat = p2/p1
3654:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

3656:       xnum = 1.0d+0 + (grat * prat)
3657:       xdenom = grat + prat
3658: 
3659:       r2 = r1 * (xnum/xdenom)
3660: 


3663:       if (debug) then
3664:         write(*,*)
3665:         write(*,*) 'r1  = ',r1
3666:         write(*,*) 'r2  = ',r2
3667:       endif

3669:       return
3670:       end
3671:       subroutine  secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err, rrg, rlg,rurg, rulg, erg, elg)
3672: !23456789012345678901234567890123456789012345678901234567890123456789012
3673: !
3674: !          subroutine secondq
3675: !
3676: !  this subroutine computes the second order (based on quick) left
3677: !  and right states for the godunov solver.
3678: !
3679: !23456789012345678901234567890123456789012345678901234567890123456789012


3682: !#######################################################################

3684:       implicit none

3686:       include 'ex74fcomd.h'

3688:       double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3689:       double precision rrg, rlg,rurg, rulg, erg, elg



3693:       double precision veld, ull,ul,ur,urr, ulg, urg

3695:       double precision fluxlim2


3698: !#######################################################################

3700: !
3701: !  compute the velocities
3702: !
3703:       ull = rull/rll
3704:       ul  = rul /rl
3705:       ur  = rur /rr
3706:       urr = rurr/rrr

3708: !
3709: !  compute the left state first
3710: !
3711:       veld = 1.0d+0

3713:       rlg = fluxlim2(rll,rl,rr,rrr,veld)
3714:       ulg = fluxlim2(ull,ul,ur,urr,veld)
3715:       rulg = rlg * ulg
3716:       elg = fluxlim2(ell,el,er,err,veld)
3717: !
3718: !  now compute the right state
3719: !
3720:       veld = -1.0d+0

3722:       rrg = fluxlim2(rll,rl,rr,rrr,veld)
3723:       urg = fluxlim2(ull,ul,ur,urr,veld)
3724:       rurg = rrg * urg
3725:       erg = fluxlim2(ell,el,er,err,veld)



3729:       return
3730:       end
3731:       double precision function shockp(x)

3733:       implicit none

3735:       double precision x
3736:       double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
3737:       logical debug

3739:       include 'ex74ftube.h'

3741:       debug = .false.


3744:       if (debug) then
3745:          write(*,*)
3746:          write(*,*) 'gamma = ',gamma
3747:          write(*,*) 'a1 = ',a1
3748:          write(*,*) 'a4 = ',a4
3749:          write(*,*) 'p1 = ',p1
3750:          write(*,*) 'p2 = ',x
3751:          write(*,*)
3752:       endif

3754:       xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
3755:       xdenom = sqrt  (  2.0d+0 * gamma * ( (2.0d+0*gamma) + (gamma + 1.0d+0) * ((x/p1) - 1) )  )
3756:       xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)

3758:       shockp = (x/p1)*((1.0d+0-(xnum/xdenom))**xpow) - (p4/p1)


3761:       if (debug) then
3762:          write(*,*)
3763:          write(*,*) 'xnum = ',xnum
3764:          write(*,*) 'gamma = ',gamma
3765:          write(*,*) 'a1 = ',a1
3766:          write(*,*) 'a4 = ',a4
3767:          write(*,*) 'p1 = ',p1
3768:          write(*,*) 'xdenom = ',xdenom
3769:          write(*,*) 'xpow = ',xpow
3770:          write(*,*) 'shockp = ',shockp
3771:          write(*,*) 'p2 = ',x
3772:          write(*,*)
3773:       endif

3775:       return
3776:       end
3777:       double precision function  uexact(x,t)

3779:       implicit none

3781:       double precision x,t
3782:       double precision xot, head, tail


3785:       logical debug

3787:       include 'ex74ftube.h'

3789:       debug = .false.

3791:       if (debug) then
3792:         write(*,*)
3793:         write(*,*) 't = ',t
3794:         write(*,*) 'x = ',x
3795:         write(*,*) 'a4 = ',a4
3796:         write(*,*) 'v3 = ',v3
3797:         write(*,*) 'a3 = ',a3
3798:         write(*,*)
3799:       endif

3801:       if (t .le. 0.0d+0) then
3802:         uexact = 0.0d+0
3803:       else

3805:        xot = x/t
3806:        head = -a4
3807:        tail = v3 - a3

3809:        if (xot .lt. head) then
3810:           uexact = 0.0d+0
3811:        elseif (xot .gt. sspd) then
3812:           uexact = 0.0d+0
3813:        elseif (xot .gt. tail) then
3814:           uexact = v2
3815:        else
3816:           uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
3817:        endif

3819:       endif


3822:       if (debug) then
3823:         write(*,*)
3824: !VAM    write(*,*) 'x = ',x,' t = ',t
3825:         write(*,*) 'uexact = ',uexact
3826:         write(*,*)
3827:       endif

3829:       return
3830:       end
3831:       double precision function upwind(fw, fe, vp)
3832: !23456789012345678901234567890123456789012345678901234567890123456789012
3833: !
3834: !          function upwind
3835: !
3836: !  this function computes the upwind face value
3837: !
3838: !23456789012345678901234567890123456789012345678901234567890123456789012


3841: !#######################################################################

3843:       implicit none

3845:       double precision fw, fe, vp

3847: !#######################################################################

3849:       if (vp .gt. 0.0) then
3850:          upwind = vp * fw
3851:       else
3852:          upwind = vp * fe
3853:       endif

3855:       return
3856:       end
3857:       subroutine uval2

3859:       implicit none

3861:       double precision prat, grat1, grat2, arat, xnum


3864:       logical debug

3866:       include 'ex74ftube.h'

3868:       debug = .false.

3870:       prat = p2/p1
3871:       grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
3872:       grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
3873:       arat = a1/gamma

3875:       xnum = sqrt ( grat2 / (prat + grat1) )

3877:       v2 = arat * (prat - 1.0d+0) * xnum

3879:       if (debug) then
3880:         write(*,*)
3881:         write(*,*) 'v2  = ',v2
3882:       endif

3884:       return
3885:       end
3886:       subroutine val3

3888:       implicit none

3890:       double precision prat, rpow, epow, p3t


3893:       logical debug

3895:       include 'ex74ftube.h'

3897:       debug = .false.


3900:       p3 = p2

3902:       prat = p3/p4

3904:       rpow = 1.0d+0 / gamma

3906:       r3 = r4 * ( prat ** rpow )

3908:       epow = (gamma - 1.0d+0) / gamma

3910:       e3 = e4 * ( (p3/p4) ** epow )

3912:       p3t = (gamma - 1.0d+0) * r3 * e3

3914:       a3 = sqrt(gamma*p3/r3)

3916:       if (debug) then
3917:         write(*,*)
3918:         write(*,*) 'a3 = ',a3
3919:         write(*,*) 'r3 = ',r3
3920:         write(*,*) 'e3 = ',e3
3921:         write(*,*) 'p3 = ',p3
3922:         write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
3923:         write(*,*)
3924:       endif

3926:       return
3927:       end
3928:       subroutine wval

3930:       implicit none

3932:       double precision prat, grat, xnum


3935:       logical debug

3937:       include 'ex74ftube.h'

3939:       debug = .false.

3941:       prat = p2/p1
3942:       grat = (gamma + 1.0d+0) / (2.0d+0 * gamma)

3944:       xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
3945: 
3946:       sspd = a1 * sqrt(xnum)
3947: 


3950:       if (debug) then
3951:         write(*,*)
3952:         write(*,*) 'sspd  = ',sspd
3953:       endif

3955:       return
3956:       end