Actual source code: ex33f.F

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

 40:       implicit none

 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !                    Include files
 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45: !
 46: !  The following include statements are generally used in SNES Fortran
 47: !  programs:
 48: !     petscsys.h  - base PETSc routines
 49: !     vec.h    - vectors
 50: !     mat.h    - matrices
 51: !     ksp.h    - Krylov subspace methods
 52: !     pc.h     - preconditioners
 53: !     snes.h   - SNES interface
 54: !  In addition, we need the following for use of PETSc drawing routines
 55: !     draw.h   - drawing routines
 56: !  Other include statements may be needed if using additional PETSc
 57: !  routines in a Fortran program, e.g.,
 58: !     viewer.h - viewers
 59: !     is.h     - index sets
 60: !
 61:  #include finclude/petscsys.h
 62:  #include finclude/petscvec.h
 63:  #include finclude/petscis.h
 64:  #include finclude/petscdraw.h
 65:  #include finclude/petscmat.h
 66:  #include finclude/petscksp.h
 67:  #include finclude/petscpc.h
 68:  #include finclude/petscsnes.h
 69:  #include finclude/petscviewer.h

 71: #include "comd.h"
 72: #include "tube.h"

 74: !
 75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76: !                   Variable declarations
 77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 78: !
 79: !  Variables:
 80: !     snes        - nonlinear solver
 81: !     x,r         - solution, residual vectors
 82: !     J           - Jacobian matrix
 83: !     its         - iterations for convergence
 84: !     matrix_free - flag - 1 indicates matrix-free version
 85: !     dt          - time step size
 86: !     draw        - drawing context
 87: !
 88:       PetscFortranAddr   ctx(6)
 89:       integer            nx,ny
 90:       SNES               snes
 91:       KSP                ksp
 92:       PC                 pc
 93:       Vec                x,r
 94:       PetscViewer        view0,view1,view2,
 95:      1                   view3, view4
 96:       Mat                Psemi
 97:       integer            matrix_free, flg, N, ierr, ngraph
 98:       integer            nstep, ndt, size, rank, i
 99:       integer            its, lits, totits, totlits
100:       integer            ndb, npc, ndamp, nwilson, ndtcon
101:       double precision   plotim
102: !      logical            pcnew

104:       double precision krtol,katol,kdtol
105:       double precision natol,nrtol,nstol
106:       integer  kmit,nmit,nmf


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

112:       external FormFunction, FormInitialGuess,FormDt,
113:      &         PCRadSetUp, PCRadApply, FormGraph,FormDampit
114:       double precision eos

116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: !  Initialize program
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

122: c
123: c  start PETSc
124: c
125:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
126:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
127:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

129:       if (size .ne. 1) then
130:          if (rank .eq. 0) then
131:             write(6,*) 'This is a uniprocessor example only!'
132:          endif
133:          SETERRQ(1,' ',ierr)
134:       endif

136: !  Initialize problem parameters

138:       debug       = .false.
139:       dampit      = .false.
140:       wilson      = .true.
141:       dtcon       = .true.
142:       pcnew       = .true.
143:       dtmax       = 1.0d+2
144:       dtmin       = 1.00d-12
145:       dt          = 1.0d-2
146:       mx          = 100
147:       nstep       = 50
148:       matrix_free = 1
149:       probnum     = 1
150:       gorder      = 1

152:       tfinal      = 1.0d+0
153:       tplot       = 0.2d+0
154:       dtgrow      = 1.05d+0
155:       tcscal      = 0.5d+0
156:       hcscal      = 0.5d+0

158:       ihod = 3
159:       ientro = 1
160:       theta = 1.00d+0
161:       pi = 3.14159d+0

163:       zero = 0.0
164:       ngraph = 10

166:       ndb = 0
167:       npc = 1

169:       damfac = 0.9d+0

171:       gamma = 1.25d+0
172:       csubv = 1.0d+0 / (gamma - 1.0d+0)

174:       v1 = 0.0d+0
175:       v4 = 0.0d+0

177:       e1 = 1.0d+0
178:       e4 = 1.0d+0

180:       r1 = 1.0d+0
181:       r4 = 2.0d+0

183:       ru1 = r1 * v1
184:       ru4 = r4 * v4

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

189:       p1 = eos(r1,ru1,et1)
190:       p4 = eos(r4,ru4,et4)

192:       a1 = sqrt(gamma*p1/r1)
193:       a4 = sqrt(gamma*p4/r4)

195:       erg0   = 1.0d+2
196:       kappa0 = 1.0d+0
197:       kappaa = -2.0d+0
198:       kappab = 13.0d+0 / 2.0d+0
199: CVAM  kappab = 2.5d+0

201: c
202: c  load the command line options
203: c
204:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
205:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
206:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nstep',nstep,flg,
207:      &                        ierr)
208:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-debug',ndb,flg,
209:      &                        ierr)
210:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-pcnew',npc,flg,
211:      &                        ierr)
212:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ihod',ihod,flg,
213:      &                        ierr)
214:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ientro',ientro,flg,
215:      &                        ierr)
216:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',
217:      &                                            theta,flg,ierr)
218:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ngraph',ngraph,flg,
219:      &                        ierr)
220:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
221:      &                            '-damfac',damfac,flg,ierr)
222:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dampit',ndamp,flg,
223:      &                        ierr)
224:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
225:      &                                   '-wilson',nwilson,flg,ierr)
226:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-gorder',gorder,flg,
227:      &                        ierr)
228:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
229:      &                            '-probnum',probnum,flg,ierr)
230:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
231:      &                            '-kappa0',kappa0,flg,ierr)
232:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
233:      &                            '-erg0',erg0,flg,ierr)
234:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dtcon',ndtcon,flg,
235:      &                        ierr)
236:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
237:      &                            '-tfinal',tfinal,flg,ierr)
238:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
239:      &                            '-tplot',tplot,flg,ierr)
240:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
241:      &                            '-dtgrow',dtgrow,flg,ierr)
242:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
243:      &                            '-tcscal',tcscal,flg,ierr)
244:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
245:      &                            '-hcscal',hcscal,flg,ierr)
246:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
247:      &                            '-dtmax',dtmax,flg,ierr)

249:       if (ndamp .eq. 1) then
250:          dampit = .true.
251:       endif

253:       if (nwilson .eq. 0) then
254:          wilson = .false.
255:       endif

257:       if (ndb .eq. 1) then
258:          debug = .true.
259:       endif

261:       if (npc .eq. 0) then
262:          pcnew = .false.
263:       endif

265:       if (ndtcon .eq. 0) then
266:          dtcon = .false.
267:       endif

269: CVAM  if (dt .ge. dtmax .or. dt .le. dtmin) then
270: CVAM     if (rank .eq. 0) write(6,*) 'DT is out of range'
271: CVAM     SETERRA(1,0,' ')
272: CVAM  endif

274:       N       = mx*neq

276:       ctx(5) = mx
277:       ctx(6) = N

279:       if (debug) then
280:         write(*,*) 'mx = ',mx
281:       endif



285: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: !  Create nonlinear solver context
287: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

289:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

291: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
292: !  Create vector data structures; set function evaluation routine
293: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

295:       CALL MatCreate(PETSC_COMM_WORLD, ctx(2), ierr)
296:       CALL MatSetSizes(ctx(2),PETSC_DECIDE,PETSC_DECIDE,mx,mx,ierr)
297:       CALL MatSetFromOptions(ctx(2),ierr)
298:       if (debug) then
299:         call MatGetSize(ctx(2),nx,ny,ierr)
300:         write(*,*) 'number of rows = ',nx,' number of col = ',ny
301:       endif
302: c
303: c  full size vectors
304: c
305:       CALL VecCreate(PETSC_COMM_WORLD,x,ierr)
306:       CALL VecSetSizes(x,PETSC_DECIDE,N,ierr)
307:       CALL VecSetType(x,VECMPI,ierr)
308:       call VecSetFromOptions(x,ierr)
309:       call VecDuplicate(x,r,ierr)
310:       call VecDuplicate(x,ctx(4),ierr)
311: c
312: c set grid
313: c
314:       dx = 2.0d+0/dfloat(mx)
315:       xl0 = -1.0d+0 -(0.5d+0 * dx)

317:       if (debug) then
318:         write(*,*) 'dx = ',dx
319:       endif
320: 

322: !  Set function evaluation routine and vector.  Whenever the nonlinear
323: !  solver needs to evaluate the nonlinear function, it will call this
324: !  routine.
325: !   - Note that the final routine argument is the user-defined
326: !     context that provides application-specific data for the
327: !     function evaluation routine.

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

331: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332: !  Customize nonlinear solver; set runtime options
333: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

338:       call SNESSetFromOptions(snes,ierr)
339: c
340: c  set the line search function to damp the newton update.
341: c
342: !     if (dampit) then
343: !       call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
344: !     endif
345: c
346: c get the linear solver info from the nonlinear solver
347: c

349:       call SNESGetKSP(snes,ksp,ierr)
350:       call KSPGetPC(ksp,pc,ierr)
351: !      call KSPGetKSP(ksp,ksp1,ierr)
352: CVAM  call KSPSetType(ksp,KSPPREONLY,ierr)
353:       call KSPSetType(ksp,KSPGMRES,ierr)

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

358:       write(*,*) 
359:       write(*,*) 
360:       write(*,*) 'Linear solver'
361:       write(*,*)
362:       write(*,*) 'rtol = ',krtol
363:       write(*,*) 'atol = ',katol
364:       write(*,*) 'dtol = ',kdtol
365:       write(*,*) 'maxits = ',kmit
366:       write(*,*)
367:       write(*,*)
368:       write(*,*) 'Nonlinear solver'
369:       write(*,*)
370:       write(*,*) 'rtol = ',nrtol
371:       write(*,*) 'atol = ',natol
372:       write(*,*) 'stol = ',nstol
373:       write(*,*) 'maxits = ',nmit
374:       write(*,*) 'max func = ',nmf
375:       write(*,*)
376:       write(*,*)

378: c
379: c  Build shell based preconditioner if flag set
380: c
381:       if (pcnew) then
382:         call PCSetType(pc,PCSHELL,ierr)
383:         call PCShellSetContext(pc,ctx,ierr)
384:         call PCShellSetSetUp(pc,PCRadSetUp,ierr)
385:         call PCShellSetApply(pc,PCRadApply,ierr)
386:       endif

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

390: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391: !  Evaluate initial guess; then solve nonlinear system.
392: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
393: c
394: c  initial counters
395: c
396:       time = 0.0d+0
397:       plotim = 0.0d+0
398:       totits = 0
399:       totlits = 0

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

406:       call FormInitialGuess(x,ierr)
407: c
408: c  open a window to plot results
409: c
410:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
411:      &                    'density',0,0,300,300,view0,ierr)
412:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
413:      &                    'velocity',320,0,300,300,view1,ierr)
414:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
415:      &                    'total energy',640,0,300,300,view2,ierr)
416:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
417:      &                    'temperature',0,360,300,300,view3,ierr)
418:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
419:      &                    'pressure',320,360,300,300,view4,ierr)
420: c
421: c  graph initial conditions
422: c
423:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)
424: c
425: c  copy x into xold
426: c
427:       call VecCopy(x,ctx(4),ierr)
428:       call FormDt(snes,x,ctx,ierr)
429: c
430: c################################
431: c
432: c  TIME STEP LOOP BEGIN
433: c
434: c################################
435: c
436:       ndt = 0

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

440:         if (debug) then
441:           write(*,*)
442:           write(*,*) 'start of time loop'
443:           write(*,*)
444:           write(*,*) 'ndt = ',ndt
445:           write(*,*) 'nstep = ',nstep
446:           write(*,*) 'time = ',time
447:           write(*,*) 'tfinal = ',tfinal
448:           write(*,*)
449:         endif

451:         ndt = ndt + 1
452: c
453: c  increment time
454: c
455:         time = time + dt
456:         plotim = plotim + dt
457: c
458: c  call the nonlinear solver
459: c
460:         call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr) 
461:         CALL SNESGetIterationNumber(snes,its,ierr)
462: c
463: c  get the number of linear iterations used by the nonlinear solver
464: c
465:         call SNESGetNumberLinearIterations(snes,lits,ierr)

467:         if (debug) then
468:            write(*,*) 'in radhyd ',ndt,'x'
469:            call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
470:         endif
471: c
472: c  increment the counters
473: c
474:         totits = totits + its
475:         totlits = totlits + lits
476: c
477: c  Compute new time step
478: c
479:           call FormDt(snes,x,ctx,ierr)
480: c
481: c  Draw contour plot of solution
482: c
483:         if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
484: 
485:            plotim = plotim - tplot


488:         if (rank .eq. 0) then
489:            write(6,100) totits,totlits,ndt,dt,time
490:         endif
491:   100   format('Newt = ',i7,' lin =',i7,' step =',i7,
492:      &         ' dt = ',e8.3,' time = ',e10.4)
493: c
494: c  graph state conditions
495: c
496:           call FormGraph(x,view0,view1,view2,view3,view4,ierr)

498:         endif
499: c
500: c copy x into xold
501: c
502:         call VecCopy(x,ctx(4),ierr)


505:         goto 10

507:       endif
508: c
509: c################################
510: c
511: c  TIME STEP LOOP END
512: c
513: c################################
514: c

516: c
517: c  graph final conditions
518: c
519:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)


522:       write(*,*)
523:       write(*,*)
524:       write(*,*) 'total Newton iterations = ',totits
525:       write(*,*) 'total linear iterations = ',totlits
526:       write(*,*) 'Average Newton per time step = ',
527:      &                       dble(totits)/dble(ndt)
528:       write(*,*) 'Average Krylov per Newton = ',
529:      &                    dble(totlits)/dble(totits)
530:       write(*,*)
531:       write(*,*)

533: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
534: !  Free work space.  All PETSc objects should be destroyed when they
535: !  are no longer needed.
536: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


539:       call MatDestroy(ctx(2),ierr)
540:       call VecDestroy(x,ierr)
541:       call VecDestroy(ctx(4),ierr)
542:       call VecDestroy(r,ierr)
543:       call SNESDestroy(snes,ierr)
544:       call PetscViewerDestroy(view0,ierr)
545:       call PetscViewerDestroy(view1,ierr)
546:       call PetscViewerDestroy(view2,ierr)
547:       call PetscViewerDestroy(view3,ierr)
548:       call PetscViewerDestroy(view4,ierr)

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

552:       call PetscFinalize(ierr)

554:       close(87)

556:       stop
557:       end
558:       subroutine ApplicationDampit(x,deltx,w,ierr)
559: ! ---------------------------------------------------------------------
560: !
561: !  ApplicationDampit - Damps the newton update, called by
562: !  the higher level routine FormDampit().
563: !
564: !  Input Parameter:
565: !  x    - current iterate
566: !  deltx   - update
567: !
568: !  Output Parameters:
569: !  x    - new iterate
570: !  ierr - error code
571: !
572: !  Notes:
573: !  This routine only damps the density.  May want to add energy
574: !  in the future
575: !

577:       implicit none

579: !  Common blocks:
580: #include "comd.h"

582: !  Input/output variables:
583:       PetscScalar x(mx*neq),deltx(mx*neq),
584:      1            w(mx*neq)
585:       integer  ierr

587: !  Local variables:
588:       double precision facmin, fac, newx, xmin, se, dse
589:       double precision u,en,rn,run
590:       integer  i, jr, jru, je

592:       0

594:       if (debug) then
595:         write(*,*) 'begin damping'
596:         do i = 1,mx*neq
597:           write(*,*)'x(',i,') = ',x(i)
598:         enddo
599:         write(*,*)
600:         do i = 1,mx*neq
601:           write(*,*)'deltx(',i,') = ',deltx(i)
602:         enddo
603:       endif

605:       facmin = 1.0d+0
606: c
607: c  set the scale factor
608: c
609:       do i=1,mx
610: c
611: c  set pointers
612: c
613:         jr  = (neq*i) - 2
614:         jru = (neq*i) - 1
615:         je  = (neq*i)
616: c
617: c  make sure dencity stayes positive
618: c
619:         newx = x(jr) - deltx(jr)
620:         xmin = damfac * x(jr)

622:         if (newx .lt. xmin) then
623:           fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
624:           if (fac .lt. facmin) then
625:             if (debug) then
626:               write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
627:             endif
628:             facmin = fac
629:           endif
630:         endif
631: c
632: c  make sure Total energy stayes positive
633: c
634:         newx = x(je) - deltx(je)
635:         xmin = damfac * x(je)

637:         if (newx .lt. xmin) then
638:           fac = (1.0d+0 - damfac)*x(je)/deltx(je)
639:           if (fac .lt. facmin) then
640:             if (debug) then
641:               write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
642:             endif
643:             facmin = fac
644:           endif
645:         endif
646: c
647: c  make sure specific internal  energy stayes positive
648: c
649: 
650:         u = x(jru)/x(jr)
651:         se = (x(je)/x(jr)) - (0.5d+0 * u * u)

653:         en = x(je) - deltx(je)
654:         rn = x(jr) - deltx(jr)
655:         run = x(jru) - deltx(jru)

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


660:         newx = se - dse
661:         xmin = damfac * se

663:         if (newx .lt. xmin) then
664:           fac = (1.0d+0 - damfac) * se / dse
665:           if (fac .lt. facmin) then
666:             if (debug) then
667:               write(*,*) 'se',i, damfac,facmin,fac,se,dse
668:             endif
669:             facmin = fac
670:           endif
671:         endif

673:       enddo
674: c
675: c write out warning
676: c
677:       if (facmin .lt. 1.0d+0) then
678:         write(*,*) 'facmin = ',facmin, damfac,time
679: c
680: c  scale the vector
681: c
682:         do i=1,neq*mx
683:           w(i) = x(i) - (facmin * deltx(i))
684:         enddo
685:       else
686:         do i=1,neq*mx
687:           w(i) = x(i) -  deltx(i)
688:         enddo
689:       endif

691:       if (debug) then
692:         write(*,*) 'end damping'
693:         do i = 1,mx*neq
694:            write(*,*) 'w(',i,') = ',w(i)
695:         enddo
696:       endif

698:       return
699:       end
700:       subroutine ApplicationDt(x,xold,ierr)
701: ! ---------------------------------------------------------------------
702: !
703: !  ApplicationDt - compute CFL numbers. Called by
704: !  the higher level routine FormDt().
705: !
706: !  Input Parameter:
707: !  x    - local vector data
708: !
709: !  Output Parameters:
710: !  ierr - error code
711: !
712: !  Notes:
713: !  This routine uses standard Fortran-style computations over a 2-dim
714: array.
715: !

717:       implicit none

719: !  Common blocks:
720: #include "comd.h"
721: #include "tube.h"

723: !  Input/output variables:
724:       PetscScalar   x(mx*neq), xold(mx*neq)
725:       integer  ierr

727: !  Local variables:
728:       integer  i, jr, jru, je
729: c
730: c new
731: c
732:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww,
733:      &                 rhouee, rhoue, rhoup, rhouw, rhouww,
734:      &                 ergee,  erge,  ergp,  ergw,  ergww,
735:      &                         vele,  velp,  velw
736:       double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
737:       double precision tcflg, hcflg, dtt, dth
738:       double precision te, tp, tw
739:       double precision ue, up, uw
740:       double precision see, sep, sew
741: c
742: c old
743: c
744:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
745:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
746:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
747:       double precision veloe,  velop,  velow
748:       double precision uop, seop, top
749:       double precision dtold, dttype
750: c
751: c  functions
752: c
753:       double precision eos

755:       dtold = dt

757:       0

759:       if (debug) then
760:         write(*,*) 'in start dt'
761:         do i = 1,mx*neq
762:           write(*,*)'x(',i,') = ',x(i)
763:         enddo
764:         write(*,*) 'tfinal = ',tfinal
765:         write(*,*) 'time = ',time
766:         write(*,*) 'dt = ',dt
767:         write(*,*) 'dtmax = ',dtmax
768:       endif

770:       sndp = -1.0d+20
771:       vradn = 0.0d+0
772:       vradd = 0.0d+0

774: c
775: c################################
776: c
777: c  loop over all cells begin
778: c
779: c################################
780: c
781:       do i=1,mx
782: c
783: c  set pointers
784: c
785:         jr  = (neq*i) - 2
786:         jru = (neq*i) - 1
787:         je  = (neq*i)
788: c
789: c
790: c  set scalars
791: c
792:         call Setpbc(i,x,
793:      &             rhoee,  rhoe,  rhop,  rhow,  rhoww,
794:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
795:      &             ergee,  erge,  ergp,  ergw,  ergww,
796:      &                     vele,  velp,  velw)
797: c
798: c compute temperatures
799: c
800:         uw = rhouw / rhow
801:         up = rhoup / rhop
802:         ue = rhoue / rhoe

804:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
805:         sep = (ergp/rhop) - (0.5d+0 * up * up)
806:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

808:         te  = see / csubv
809:         tp  = sep / csubv
810:         tw  = sew / csubv
811: c
812: c compute old temperature
813: c

815:         call Setpbc(i,xold,
816:      &             rhooee,  rhooe,  rhoop,  rhoow,  rhooww,
817:      &             rhouoee, rhouoe, rhouop, rhouow, rhouoww,
818:      &             ergoee,  ergoe,  ergop,  ergow,  ergoww,
819:      &                      veloe,  velop,  velow)

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

826:         uop = rhouop / rhoop

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

830:         top  = seop / csubv
831: c
832: c  compute thermal cfl
833: c
834:         vradn = vradn + (abs(tp - top)/dt)
835:         vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
836: c
837: c  compute hydro cfl
838: c

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

843:       enddo
844: c
845: c################################
846: c
847: c  loop over all cells end
848: c
849: c################################
850: c

852:       vrad = vradn / vradd

854:       tcfl = (vrad * dt) / dx
855:       hcfl = (sndp * dt) / dx

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

860:       dth = hcscal * dx / sndp

862:       if (.not. dtcon) then
863:         dt = min (dth,dtt,dt*dtgrow)
864:         if (dt .lt. dtmin) then
865:            dt = dtmin
866:         endif
867:         if (dt .gt. dtmax) then
868:            dt = dtmax
869:         endif
870:         if ( (time + dt) .gt. tfinal) then
871:            dt = tfinal - time
872:         endif

874:         if (dt .eq. dth) then
875:            dttype = 1.0d+0
876:         elseif (dt .eq. dtt) then
877:            dttype = 2.0d+0
878:         elseif (dt .eq. dtold*dtgrow) then
879:            dttype = 3.0d+0
880:         elseif (dt .eq. dtmax) then
881:            dttype = 4.0d+0
882:         elseif (dt .eq. dtmin) then
883:            dttype = 5.0d+0
884:         elseif (dt .eq. tfinal - time) then
885:            dttype = 6.0
886:         else
887:            dttype = -1.0d+0
888:         endif

890:       endif
891: 
892: 
893:       write(87,1000) time,dt,dth/hcscal,dtt/tcscal
894:       write(88,1000) time,dttype

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

898:       if (debug) then
899:         write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
900:         write(*,*) 'dtt = ',dtt,' dth = ',dth
901:         write(*,*) 'tfinal = ',tfinal
902:         write(*,*) 'time = ',time
903:         write(*,*) 'dt = ',dt
904:         write(*,*) 'dtmax = ',dtmax
905:         write(*,*)
906:       endif

908:       return
909:       end
910:       subroutine ApplicationExact(x,ierr)
911: ! ---------------------------------------------------------------------
912: !
913: !  ApplicationExact - Computes exact solution, called by
914: !  the higher level routine FormExact().
915: !
916: !  Input Parameter:
917: !  x - local vector data
918: !
919: !  Output Parameters:
920: !  x -    initial conditions
921: !  ierr - error code
922: !
923: !  Notes:
924: !  This routine uses standard Fortran-style computations over a 1-dim
925: array.
926: !

928:       implicit none

930: !  Common blocks:

932: #include "comd.h"

934: !  Input/output variables:
935:       PetscScalar  x(mx)
936:       integer ierr

938: !  Local variables:
939:       integer  i
940:       double precision xloc
941:       PetscScalar rexact


944: !  Set parameters

946:       0

948:       do i = 1,mx

950:         xloc = xl0 + (dble(i) * dx)
951:         x(i) = rexact(xloc,time)

953:       enddo

955:       return
956:       end
957:       subroutine ApplicationFunction(x,f,xold,ierr)
958: ! ---------------------------------------------------------------------
959: !
960: !  ApplicationFunction - Computes nonlinear function, called by
961: !  the higher level routine FormFunction().
962: !
963: !  Input Parameter:
964: !  x    - local vector data
965: !
966: !  Output Parameters:
967: !  f    - local vector data, f(x)
968: !  ierr - error code
969: !
970: !  Notes:
971: !  This routine uses standard Fortran-style computations over a 2-dim
972: array.
973: !

975:       implicit none

977: !  Common blocks:
978: #include "comd.h"

980: !  Input/output variables:
981:       PetscScalar   x(mx*neq),f(mx*neq),
982:      1              xold(mx*neq)
983:       integer       ierr

985: !  Local variables:
986:       integer  i, jr, jru, je
987:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww,
988:      &                 rhouee, rhoue, rhoup, rhouw, rhouww,
989:      &                 ergee,  erge,  ergp,  ergw,  ergww,
990:      &                         vele,  velp,  velw

992:       double precision cont, energy, mom

994:       0

996:       if (debug) then
997:         write(*,*) 'in function'
998:         do i = 1,mx*neq
999:           write(*,*)'x(',i,') = ',x(i)
1000:         enddo
1001:       endif
1002: c
1003: c################################
1004: c
1005: c  loop over all cells begin
1006: c
1007: c################################
1008: c
1009:       do i=1,mx
1010: c
1011: c  set pointers
1012: c
1013:       jr  = (neq*i) - 2
1014:       jru = (neq*i) - 1
1015:       je  = (neq*i)
1016: c
1017: c
1018: c  set scalars
1019: c
1020:       call Setpbc(i,x,
1021:      &             rhoee,  rhoe,  rhop,  rhow,  rhoww,
1022:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
1023:      &             ergee,  erge,  ergp,  ergw,  ergww,
1024:      &                     vele,  velp,  velw)
1025: c
1026: c  compute functions
1027: c

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


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


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

1046:        if (debug) then
1047:          write(*,*)
1048:          write(*,*) i,jr,jru,je,'res,r,ru,e'
1049:          write(*,*) f(jr),f(jru),f(je)
1050:          write(*,*)
1051:        endif

1053:       enddo
1054: c
1055: c################################
1056: c
1057: c  loop over all cells end
1058: c
1059: c################################
1060: c

1062:       if (debug) then
1063:         write(*,*) 'in function'
1064:         do i = 1,mx*neq
1065:            write(*,*) 'f(',i,') = ',f(i)
1066:         enddo
1067:       endif

1069:       return
1070:       end
1071:       subroutine ApplicationInitialGuess(x,ierr)
1072: ! ---------------------------------------------------------------------
1073: !
1074: !  ApplicationInitialGuess - Computes initial approximation, called by
1075: !  the higher level routine FormInitialGuess().
1076: !
1077: !  Input Parameter:
1078: !  x - local vector data
1079: !
1080: !  Output Parameters:
1081: !  x -    initial conditions
1082: !  ierr - error code
1083: !
1084: !  Notes:
1085: !  This routine uses standard Fortran-style computations over a 1-dim
1086: array.
1087: !

1089:       implicit none

1091: !  Common blocks:

1093: #include "comd.h"
1094: #include "tube.h"

1096: !  Input/output variables:
1097:       PetscScalar  x(mx*neq)
1098:       integer ierr

1100: !  Local variables:
1101:       integer  i, j, jr, jru, je
1102:       double precision xloc, re, ee, ve
1103:       double precision wid, efloor
1104:       PetscScalar uexact, rexact, eexact


1107: CVAM  efloor = max(1.0d-1,1.0d-3 * erg0)
1108:       efloor = 1.0d-1
1109: CVAM  wid = max(1.0d-2,dx)
1110:       wid = 1.0d-2

1112: !  Set parameters

1114:       0

1116:       do i = 1,mx

1118:         jr  = (neq*i) - 2
1119:         jru = (neq*i) - 1
1120:         je  = (neq*i)

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

1124:         if (probnum .eq. 1) then
1125:            re = rexact(xloc,zero)
1126:            ve = uexact(xloc,zero)
1127:            ee = eexact(xloc,zero)
1128:         else
1129:            re = 1.0d+0
1130:            ve = 0.0d+0
1131:            ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1132:         endif

1134:         x(jr)  = re
1135:         x(jru) = re * ve
1136:         x(je)  = re * ( (0.5d+0 * ve * ve) + ee )

1138:         if (debug) then
1139:            write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1140:  100       format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1141:         endif

1143:       enddo

1145:       call exact0
1146:       call eval2
1147:       call rval2
1148:       call wval
1149:       call uval2
1150:       v3 = v2
1151:       call val3

1153:       a1 = sqrt(gamma*p1/r1)
1154:       a2 = sqrt(gamma*p2/r2)
1155:       a3 = sqrt(gamma*p3/r3)
1156:       a4 = sqrt(gamma*p4/r4)

1158:       write(*,1000) r1,r2,r3,r4
1159:       write(*,2000) p1,p2,p3,p4
1160:       write(*,3000) e1,e2,e3,e4
1161:       write(*,4000) a1,a2,a3,a4
1162:       write(*,*)

1164:  1000 format ('rhos      ',4(f12.6))
1165:  2000 format ('pressures ',4(f12.6))
1166:  3000 format ('energies  ',4(f12.6))
1167:  4000 format ('sound     ',4(f12.6))


1170:       return
1171:       end
1172:       subroutine ApplicationXmgr(x,ivar,ierr)
1173: ! ---------------------------------------------------------------------
1174: !
1175: !  ApplicationXmgr - Sets the Xmgr output called from
1176: !  the higher level routine FormXmgr().
1177: !
1178: !  Input Parameter:
1179: !  x - local vector data
1180: !
1181: !  Output Parameters:
1182: !  x -    initial conditions
1183: !  ierr - error code
1184: !
1185: !  Notes:
1186: !  This routine uses standard Fortran-style computations over a 1-dim
1187: array.
1188: !

1190:       implicit none

1192: !  Common blocks:

1194: #include "comd.h"

1196: !  Input/output variables:
1197:       PetscScalar  x(mx)
1198:       integer ivar,ierr

1200: !  Local variables:
1201:       integer  i
1202:       double precision xloc, sum
1203:       PetscScalar rexact
1204:       integer iplotnum(5)
1205:       save iplotnum
1206:       character*8 grfile


1209:       data iplotnum / -1,-1,-1,-1,-1 /



1213: !  Set parameters

1215:       iplotnum(ivar) = iplotnum(ivar) + 1
1216:       0

1218:       if (ivar .eq. 1) then
1219:          write(grfile,4000) iplotnum(ivar)
1220:  4000    format('Xmgrr',i3.3)
1221:       elseif (ivar .eq. 2) then
1222:          write(grfile,5000) iplotnum(ivar)
1223:  5000    format('Xmgru',i3.3)
1224:       elseif (ivar .eq. 3) then
1225:          write(grfile,6000) iplotnum(ivar)
1226:  6000    format('Xmgre',i3.3)
1227:       elseif (ivar .eq. 4) then
1228:          write(grfile,7000) iplotnum(ivar)
1229:  7000    format('Xmgrt',i3.3)
1230:       else
1231:          write(grfile,8000) iplotnum(ivar)
1232:  8000    format('Xmgrp',i3.3)
1233:       endif

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

1237:       do i = 1,mx

1239:         xloc = xl0 + (dble(i) * dx)
1240:         if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1241:           write(44,1000) xloc, x(i), rexact(xloc,time)
1242:         else
1243:           write(44,1000) xloc, x(i)
1244:         endif

1246:       enddo

1248:  1000 format(3(e18.12,2x))
1249:       close(44)

1251:       if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1252:         sum = 0.0d+0
1253:         do i = 1,mx
1254:            xloc = xl0 + (dble(i) * dx)
1255:            sum = sum + (x(i) - rexact(xloc,time)) ** 2
1256:         enddo
1257:         sum = sqrt(sum)

1259:         write(*,*)
1260:         write(*,*)  'l2norm of the density error is',sum
1261:         write(*,*)
1262:       endif


1265:       return
1266:       end
1267:       subroutine FormDampit(snes,ctx,x,f,g,y,w,
1268:      &                       fnorm,ynorm,gnorm,flag,ierr)
1269: ! ---------------------------------------------------------------------
1270: !
1271: !  FormDampit - damps the Newton update
1272: !
1273: !  Input Parameters:
1274: !  snes  - the SNES context
1275: !  x     - current iterate
1276: !  f     - residual evaluated at x
1277: !  y     - search direction (containes new iterate on output)
1278: !  w     - work vector
1279: !  fnorm - 2-norm of f
1280: !
1281: !  In this example the application context is a Fortran integer array:
1282: !      ctx(1) = shell preconditioner pressure matrix contex
1283: !      ctx(2) = semi implicit pressure matrix
1284: !      ctx(4) = xold  - old time values need for time advancement
1285: !      ctx(5) = mx    - number of control volumes
1286: !      ctx(6) = N     - total number of unknowns
1287: !
1288: !  Output Parameter:
1289: !  g     - residual evaluated at new iterate y
1290: !  y     - new iterate (contains search direction on input
1291: !  gnorm - 2-norm of g
1292: !  ynorm - 2-norm of search length
1293: !  flag  - set to 0 if the line search succeeds; -1 on failure
1294: !
1295: !  Notes:
1296: !  This routine serves as a wrapper for the lower-level routine
1297: !  "ApplicationDampit", where the actual computations are
1298: !  done using the standard Fortran style of treating the local
1299: !  vector data as a multidimensional array over the local mesh.
1300: !  This routine merely accesses the local vector data via
1301: !  VecGetArray() and VecRestoreArray().
1302: !
1303:       implicit none

1305:  #include finclude/petscsys.h
1306:  #include finclude/petscvec.h
1307:  #include finclude/petscsnes.h

1309: !  Input/output variables:
1310:       SNES             snes
1311:       Vec              x, f, g, y, w
1312:       PetscFortranAddr ctx(*)
1313:       PetscScalar           fnorm, ynorm, gnorm
1314:       integer          ierr, flag

1316: !  Common blocks:

1318: #include "comd.h"

1320: !  Local variables:

1322: !  Declarations for use with local arrays:
1323:       PetscScalar lx_v(0:1),ly_v(0:1),
1324:      1            lw_v(0:1)
1325:       PetscOffset lx_i,ly_i,lw_i

1327: c
1328: c  set ynorm
1329: c
1330:       call VecNorm(y,NORM_2,ynorm,ierr)
1331: c
1332: c  copy x to w
1333: c
1334:       call VecCopy(x,w,ierr)
1335: c
1336: c get pointers to x, y, w
1337: c

1339:       call VecGetArray(x,lx_v,lx_i,ierr)
1340:       call VecGetArray(y,ly_v,ly_i,ierr)
1341:       call VecGetArray(w,lw_v,lw_i,ierr)
1342: c
1343: c  Compute Damping 
1344: c
1345:       call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1346: c
1347: c  Restore vectors x, y, w
1348: c
1349:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1350:       call VecRestoreArray(y,ly_v,ly_i,ierr)
1351:       call VecRestoreArray(w,lw_v,lw_i,ierr)
1352: c
1353: c  copy w to y
1354: c
1355:       call VecCopy(w,y,ierr)
1356: c
1357: c  compute new residual
1358: c
1359:       call SNESComputeFunction(snes,y,g,ierr)
1360:       call VecNorm(g,NORM_2,gnorm,ierr)
1361:       flag = 0

1363:       if (debug) then
1364:          write(*,*) 'form damp ynorm = ',ynorm
1365:          write(*,*) 'gnorm = ',gnorm
1366:       endif

1368:       return
1369:       end
1370:       subroutine FormDt(snes,x,ctx,ierr)
1371: ! ---------------------------------------------------------------------
1372: !
1373: !  FormDt - Compute CFL numbers
1374: !
1375: !  Input Parameters:
1376: !  snes  - the SNES context
1377: !  x     - input vector
1378: !
1379: !  In this example the application context is a Fortran integer array:
1380: !      ctx(1) = shell preconditioner pressure matrix contex
1381: !      ctx(2) = semi implicit pressure matrix
1382: !      ctx(4) = xold  - old time values need for time advancement
1383: !      ctx(5) = mx    - number of control volumes
1384: !      ctx(6) = N     - total number of unknowns
1385: !
1386: !
1387: !  Notes:
1388: !  This routine serves as a wrapper for the lower-level routine
1389: !  "ApplicationDt", where the actual computations are
1390: !  done using the standard Fortran style of treating the local
1391: !  vector data as a multidimensional array over the local mesh.
1392: !  This routine merely accesses the local vector data via
1393: !  VecGetArray() and VecRestoreArray().
1394: !
1395:       implicit none

1397:  #include finclude/petscsys.h
1398:  #include finclude/petscvec.h
1399:  #include finclude/petscsnes.h

1401: !  Input/output variables:
1402:       SNES             snes
1403:       Vec              x
1404:       PetscFortranAddr ctx(*)
1405:       integer          ierr

1407: !  Common blocks:

1409: #include "comd.h"

1411: !  Local variables:

1413: !  Declarations for use with local arrays:
1414:       PetscScalar      lx_v(0:1)
1415:       PetscOffset lx_i
1416:       PetscScalar      lxold_v(0:1)
1417:       PetscOffset lxold_i

1419: c
1420: c get pointers to x, xold
1421: c

1423:       call VecGetArray(x,lx_v,lx_i,ierr)
1424:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1425: c
1426: c  Compute function 
1427: c
1428:       call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1429: c
1430: c  Restore vectors x, xold
1431: c
1432:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1433:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)

1435:       return
1436:       end
1437:       subroutine FormExact(x,ierr)
1438: ! ---------------------------------------------------------------------
1439: !
1440: !  FormExact - Forms exact solution
1441: !
1442: !  Input Parameter:
1443: !  x - vector
1444: !
1445: !  Output Parameters:
1446: !  x - vector
1447: !  ierr - error code
1448: !
1449: !  Notes:
1450: !  This routine serves as a wrapper for the lower-level routine
1451: !  "ApplicationExact", where the actual computations are
1452: !  done using the standard Fortran style of treating the local
1453: !  vector data as a multidimensional array over the local mesh.
1454: !  This routine merely accesses the local vector data via
1455: !  VecGetArray() and VecRestoreArray().
1456: !
1457:       implicit none

1459:  #include finclude/petscsys.h
1460:  #include finclude/petscvec.h
1461:  #include finclude/petscmat.h
1462:  #include finclude/petscsnes.h

1464: !  Input/output variables:
1465:       Vec      x
1466:       integer  ierr

1468: !  Declarations for use with local arrays:
1469:       PetscScalar      lx_v(0:1)
1470:       PetscOffset lx_i

1472:       0

1474: c
1475: c  get a pointer to x
1476: c
1477:       call VecGetArray(x,lx_v,lx_i,ierr)
1478: c
1479: c  Compute initial guess 
1480: c
1481:       call ApplicationExact(lx_v(lx_i),ierr)
1482: c
1483: c  Restore vector x
1484: c
1485:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1487:       return 
1488:       end
1489:       subroutine FormFunction(snes,x,f,ctx,ierr)
1490: ! ---------------------------------------------------------------------
1491: !
1492: !  FormFunction - Evaluates nonlinear function, f(x).
1493: !
1494: !  Input Parameters:
1495: !  snes  - the SNES context
1496: !  x     - input vector
1497: !
1498: !  In this example the application context is a Fortran integer array:
1499: !      ctx(1) = shell preconditioner pressure matrix contex
1500: !      ctx(2) = semi implicit pressure matrix
1501: !      ctx(4) = xold  - old time values need for time advancement
1502: !      ctx(5) = mx    - number of control volumes
1503: !      ctx(6) = N     - total number of unknowns
1504: !
1505: !  Output Parameter:
1506: !  f     - vector with newly computed function
1507: !
1508: !  Notes:
1509: !  This routine serves as a wrapper for the lower-level routine
1510: !  "ApplicationFunction", where the actual computations are
1511: !  done using the standard Fortran style of treating the local
1512: !  vector data as a multidimensional array over the local mesh.
1513: !  This routine merely accesses the local vector data via
1514: !  VecGetArray() and VecRestoreArray().
1515: !
1516:       implicit none

1518:  #include finclude/petscsys.h
1519:  #include finclude/petscvec.h
1520:  #include finclude/petscsnes.h

1522: !  Input/output variables:
1523:       SNES             snes
1524:       Vec              x, f
1525:       PetscFortranAddr ctx(*)
1526:       integer          ierr

1528: !  Common blocks:

1530: #include "comd.h"

1532: !  Local variables:

1534: !  Declarations for use with local arrays:
1535:       PetscScalar      lx_v(0:1), lf_v(0:1)
1536:       PetscOffset lx_i, lf_i
1537:       PetscScalar      lxold_v(0:1)
1538:       PetscOffset lxold_i

1540: c
1541: c get pointers to x, f, xold
1542: c

1544:       call VecGetArray(x,lx_v,lx_i,ierr)
1545:       call VecGetArray(f,lf_v,lf_i,ierr)
1546:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1547: c
1548: c  Compute function 
1549: c
1550:       call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),
1551:      &                             lxold_v(lxold_i),ierr)
1552: c
1553: c  Restore vectors x, f, xold
1554: c
1555:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1556:       call VecRestoreArray(f,lf_v,lf_i,ierr)
1557:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1558: c
1559: c something to do with profiling
1560: c
1561:       call PetscLogFlops(11.0d0*mx,ierr)

1563:       return
1564:       end
1565:       subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1566: ! ---------------------------------------------------------------------
1567: !
1568: !  FormGraph - Forms Graphics output
1569: !
1570: !  Input Parameter:
1571: !  x - vector
1572: !  time - scalar
1573: !
1574: !  Output Parameters:
1575: !  ierr - error code
1576: !
1577: !  Notes:
1578: !  This routine serves as a wrapper for the lower-level routine
1579: !  "ApplicationXmgr", where the actual computations are
1580: !  done using the standard Fortran style of treating the local
1581: !  vector data as a multidimensional array over the local mesh.
1582: !  This routine merely accesses the local vector data via
1583: !  VecGetArray() and VecRestoreArray().
1584: !
1585:       implicit none

1587:  #include finclude/petscsys.h
1588:  #include finclude/petscvec.h
1589:  #include finclude/petscmat.h
1590:  #include finclude/petscsnes.h

1592: #include "comd.h"
1593: #include "tube.h"

1595: !  Input/output variables:
1596:       Vec      x
1597:       integer  ierr

1599: !  Declarations for use with local arrays:
1600:       IS                 rfrom,rto,
1601:      1                   rufrom, ruto, efrom, eto
1602:       Vec                rval
1603:       Vec                uval
1604:       Vec                ruval
1605:       Vec                eval
1606:       Vec                seval
1607:       Vec                pval
1608:       Vec                kval
1609:       Vec                tval
1610:       Vec                steval
1611:       VecScatter         scatter
1612:       PetscViewer        view0,view1,
1613:      1                   view2, view3, view4
1614:       double precision   minus1, l2err, gm1, csubvi


1617:       csubvi = 1.0d+0 / csubv
1618:       gm1 = gamma - 1.0d+0
1619:       0
1620:       minus1 = -1.0d+0
1621: c
1622: c  graphics vectors
1623: c
1624:       CALL VecCreate(PETSC_COMM_WORLD,rval,ierr)
1625:       CALL VecSetSizes(rval,PETSC_DECIDE,mx,ierr)
1626:       CALL VecSetType(rval,VECMPI,ierr)
1627:       call VecSetFromOptions(rval,ierr)
1628:       call VecDuplicate(rval,uval,ierr)
1629:       call VecDuplicate(rval,ruval,ierr)
1630:       call VecDuplicate(rval,eval,ierr)
1631:       call VecDuplicate(rval,seval,ierr)
1632:       call VecDuplicate(rval,pval,ierr)
1633:       call VecDuplicate(rval,kval,ierr)
1634:       call VecDuplicate(rval,tval,ierr)
1635:       call VecDuplicate(rval,steval,ierr)
1636: c
1637: c create index sets for vector scatters
1638: c
1639:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1640:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  rto,   ierr)
1641:       call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1642:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  ruto,  ierr)
1643:       call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1644:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  eto,   ierr)

1646: c
1647: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1648: c
1649: c
1650: c  load rval from x
1651: c
1652:       call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1653:       call VecScatterBegin(scatter,x,rval,INSERT_VALUES,
1654:      &                 SCATTER_FORWARD,ierr)
1655:       call VecScatterEnd(scatter,x,rval,INSERT_VALUES,
1656:      &                 SCATTER_FORWARD,ierr)
1657:       call VecScatterDestroy(scatter,ierr)
1658: c
1659: c  plot rval vector
1660: c
1661:       call VecView(rval,view0,ierr)
1662: c
1663: c  make xmgr plot of rval
1664: c
1665:       call FormXmgr(rval,1,ierr)
1666: c
1667: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1668: c
1669: c
1670: c  load eval from x
1671: c
1672:       call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1673:       call VecScatterBegin(scatter,x,eval,INSERT_VALUES,
1674:      &                 SCATTER_FORWARD,ierr)
1675:       call VecScatterEnd(scatter,x,eval,INSERT_VALUES,
1676:      &                 SCATTER_FORWARD,ierr)
1677:       call VecScatterDestroy(scatter,ierr)
1678: c
1679: c  plot eval vector
1680: c
1681:       call VecView(eval,view2,ierr)
1682: c
1683: c  make xmgr plot of eval
1684: c
1685:       call FormXmgr(eval,3,ierr)
1686: c
1687: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1688: c

1690: c
1691: c  load ruval from x
1692: c
1693:       call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1694:       call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,
1695:      &                 SCATTER_FORWARD,ierr)
1696:       call VecScatterEnd(scatter,x,ruval,INSERT_VALUES,
1697:      &                 SCATTER_FORWARD,ierr)
1698:       call VecScatterDestroy(scatter,ierr)
1699: c
1700: c  create u = ru / r
1701: c
1702:       call VecPointwiseDivide(ruval,rval,uval,ierr)
1703: c
1704: c  plot uval vector
1705: c
1706:       call VecView(uval,view1,ierr)
1707: c
1708: c  make xmgr plot of uval
1709: c
1710:       call FormXmgr(uval,2,ierr)

1712: c
1713: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1714: c

1716:       call VecPointwiseMult(kval,uval,uval,ierr)
1717:       call VecScale(kval,0.5d+0,ierr)

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

1722:       call VecCopy(seval,tval,ierr)
1723:       call VecScale(tval,csubvi,ierr)

1725: c
1726: c  plot tval vector
1727: c
1728:       call VecView(tval,view3,ierr)
1729: c
1730: c  make xmgr plot of tval
1731: c
1732:       call FormXmgr(tval,4,ierr)

1734: c
1735: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736: c

1738:       call VecPointwiseMult(rval,seval,pval,ierr)
1739:       call VecScale(pval,gm1,ierr)
1740: c
1741: c  plot pval vector
1742: c
1743:       call VecView(pval,view4,ierr)
1744: c
1745: c  make xmgr plot of pval
1746: c
1747:       call FormXmgr(pval,5,ierr)
1748: c
1749: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1750: c





1756: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1757: !  Free work space.  All PETSc objects should be destroyed when they
1758: !  are no longer needed.
1759: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1761:       call VecDestroy(rval, ierr)
1762:       call VecDestroy(uval, ierr)
1763:       call VecDestroy(ruval,ierr)
1764:       call VecDestroy(eval, ierr)
1765:       call VecDestroy(seval, ierr)
1766:       call VecDestroy(pval, ierr)
1767:       call VecDestroy(kval, ierr)
1768:       call VecDestroy(tval, ierr)
1769:       call VecDestroy(steval, ierr)

1771:       call ISDestroy(rfrom, ierr)
1772:       call ISDestroy(rto,   ierr)

1774:       call ISDestroy(rufrom,ierr)
1775:       call ISDestroy(ruto,  ierr)

1777:       call ISDestroy(efrom, ierr)
1778:       call ISDestroy(eto,   ierr)


1781:       return
1782:       end
1783:       subroutine FormInitialGuess(x,ierr)
1784: ! ---------------------------------------------------------------------
1785: !
1786: !  FormInitialGuess - Forms initial approximation.
1787: !
1788: !  Input Parameter:
1789: !  x - vector
1790: !
1791: !  Output Parameters:
1792: !  x - vector
1793: !  ierr - error code
1794: !
1795: !  Notes:
1796: !  This routine serves as a wrapper for the lower-level routine
1797: !  "ApplicationInitialGuess", where the actual computations are
1798: !  done using the standard Fortran style of treating the local
1799: !  vector data as a multidimensional array over the local mesh.
1800: !  This routine merely accesses the local vector data via
1801: !  VecGetArray() and VecRestoreArray().
1802: !
1803:       implicit none

1805:  #include finclude/petscsys.h
1806:  #include finclude/petscvec.h
1807:  #include finclude/petscmat.h
1808:  #include finclude/petscsnes.h

1810: !  Input/output variables:
1811:       Vec      x
1812:       integer  ierr

1814: !  Declarations for use with local arrays:
1815:       PetscScalar      lx_v(0:1)
1816:       PetscOffset lx_i

1818:       0

1820: c
1821: c  get a pointer to x
1822: c
1823:       call VecGetArray(x,lx_v,lx_i,ierr)
1824: c
1825: c  Compute initial guess 
1826: c
1827:       call ApplicationInitialGuess(lx_v(lx_i),ierr)
1828: c
1829: c  Restore vector x
1830: c
1831:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1833:       return 
1834:       end
1835:       subroutine FormXmgr(x,ivar,ierr)
1836: ! ---------------------------------------------------------------------
1837: !
1838: !  FormXmgr - Forms Xmgr output
1839: !
1840: !  Input Parameter:
1841: !  x - vector
1842: !
1843: !  Output Parameters:
1844: !  x - vector
1845: !  ierr - error code
1846: !
1847: !  Notes:
1848: !  This routine serves as a wrapper for the lower-level routine
1849: !  "ApplicationXmgr", where the actual computations are
1850: !  done using the standard Fortran style of treating the local
1851: !  vector data as a multidimensional array over the local mesh.
1852: !  This routine merely accesses the local vector data via
1853: !  VecGetArray() and VecRestoreArray().
1854: !
1855:       implicit none

1857:  #include finclude/petscsys.h
1858:  #include finclude/petscvec.h
1859:  #include finclude/petscmat.h
1860:  #include finclude/petscsnes.h

1862: !  Input/output variables:
1863:       Vec      x
1864:       integer  ivar,ierr

1866: !  Declarations for use with local arrays:
1867:       PetscScalar      lx_v(0:1)
1868:       PetscOffset lx_i

1870:       0

1872: c
1873: c  get a pointer to x
1874: c
1875:       call VecGetArray(x,lx_v,lx_i,ierr)
1876: c
1877: c  make the graph
1878: c
1879:       call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1880: c
1881: c  Restore vector x
1882: c
1883:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1885:       return 
1886:       end
1887:       subroutine PCRadApply(pc,x,y,ierr)
1888: ! -------------------------------------------------------------------
1889: !
1890: !   PCRadApply - This routine demonstrates the use of a
1891: !   user-provided preconditioner.
1892: !
1893: !   Input Parameters:
1894: !   pc - preconditioner object
1895: !   x - input vector
1896: !  In this example the shell preconditioner application context
1897: !  is a Fortran integer array:
1898: !      ctx(1) = shell preconditioner pressure matrix contex
1899: !      ctx(2) = semi implicit pressure matrix
1900: !      ctx(4) = xold  - old time values need for time advancement
1901: !      ctx(5) = mx    - number of control volumes
1902: !      ctx(6) = N     - total number of unknowns
1903: !
1904: !   Output Parameters:
1905: !   y - preconditioned vector
1906: !   ierr  - error code (nonzero if error has been detected)
1907: !
1908: !   Notes:
1909: !   This code implements the Jacobi preconditioner plus the
1910: !   SOR preconditioner
1911: !

1913:       implicit none

1915:  #include finclude/petscsys.h
1916:  #include finclude/petscvec.h

1918: #include "comd.h"
1919: c
1920: c  Input
1921: c
1922:       PC               pc
1923:       Vec              x, y
1924:       integer          ierr
1925: c
1926: c  Local
1927: c
1928:       PetscFortranAddr ctx(*)
1929:       IS               defrom, deto
1930:       Vec              de, rese
1931:       VecScatter       scatter
1932:       PetscScalar      lde_v(0:1),lrese_v(0:1)
1933:       PetscOffset      lde_i,     lrese_i

1935:       call PCShellGetContext(pc,ctx,ierr)

1937: c
1938: c  Identity preconditioner
1939: c
1940:       call VecCopy(x,y,ierr)
1941: c
1942: c  if kappa0 not equal to zero then precondition the radiation diffusion
1943: c
1944:       if (kappa0 .ne. 0.0d+0) then
1945: 

1947: c
1948: c  Create needed vectors
1949: c
1950:          CALL VecCreate(PETSC_COMM_WORLD,de,ierr)
1951:          CALL VecSetSizes(de,PETSC_DECIDE,mx,ierr)
1952:          CALL VecSetType(de,VECMPI,ierr)
1953:          call VecSetFromOptions(de,ierr)
1954:          call VecDuplicate(de,rese,ierr)
1955: c
1956: c  create index sets for scatters
1957: c
1958:          call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1959:          call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1960: c
1961: c  load rese from x
1962: c
1963:          call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1964:          call VecScatterBegin(scatter,x,rese,INSERT_VALUES,
1965:      &                 SCATTER_FORWARD,ierr)
1966:          call VecScatterEnd(scatter,x,rese,INSERT_VALUES,
1967:      &                 SCATTER_FORWARD,ierr)
1968:          call VecScatterDestroy(scatter,ierr)
1969: c
1970: c  apply preconditioner
1971: c
1972:       call PCApply(ctx(1),rese,de,ierr)

1974:       if (debug) then
1975:         write(*,*) 'PCRadApply dh is'
1976:         call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1977:       endif
1978: c
1979: c load de into y
1980: c
1981:       call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1982:       call VecScatterBegin(scatter,de,y,INSERT_VALUES,
1983:      &                 SCATTER_FORWARD,ierr)
1984:       call VecScatterEnd(scatter,de,y,INSERT_VALUES,
1985:      &                 SCATTER_FORWARD,ierr)
1986:       call VecScatterDestroy(scatter,ierr)

1988:       if (debug) then
1989:         write(*,*) 'PCRadApply y is'
1990:         call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1991:       endif

1993:       call VecDestroy(de,ierr)
1994:       call VecDestroy(rese,ierr)

1996:       call ISDestroy(defrom,ierr)
1997:       call ISDestroy(deto,ierr)

1999:       endif


2002:       return
2003:       end
2004:       subroutine PCRadSetUp(pc,ierr)
2005: !
2006: !   PCRadSetUp - This routine sets up a user-defined
2007: !   preconditioner context.
2008: !
2009: !   Input Parameters:
2010: !  In this example the shell preconditioner application context
2011: !  is a Fortran integer array:
2012: !      ctx(1) = shell preconditioner pressure matrix contex
2013: !      ctx(2) = semi implicit pressure matrix
2014: !      ctx(4) = xold  - old time values need for time advancement
2015: !      ctx(5) = mx    - number of control volumes
2016: !      ctx(6) = N     - total number of unknowns
2017: !
2018: !   Output Parameter:
2019: !   ierr  - error code (nonzero if error has been detected)
2020: !
2021: !   Notes:
2022: !   In this example, we define the shell preconditioner to be Jacobi
2023: !   method.  Thus, here we create a work vector for storing the reciprocal
2024: !   of the diagonal of the preconditioner matrix; this vector is then
2025: !   used within the routine PCRadApply().
2026: !

2028:       implicit none

2030:  #include finclude/petscsys.h
2031:  #include finclude/petscvec.h
2032:  #include finclude/petscmat.h

2034: #include "comd.h"
2035: c
2036: c  Input
2037: c
2038:       PC               pc
2039:       integer          ierr
2040: c
2041: c  Local
2042: c
2043:       PetscFortranAddr ctx(*)
2044:       Vec              eold
2045: 
2046:       PetscScalar      le_v(0:1)
2047:       PetscOffset le_i

2049:       call PCShellGetContext(pc,ctx,ierr)
2050: 
2051: c
2052: c  create vector
2053: c
2054:       CALL VecCreate(PETSC_COMM_WORLD,eold,ierr)
2055:       CALL VecSetSizes(eold,PETSC_DECIDE,mx,ierr)
2056:       CALL VecSetType(eold,VECMPI,ierr)
2057:       call VecSetFromOptions(eold,ierr)
2058: c
2059: c set up the matrix based on xold
2060: c
2061:       call Setmat(ctx,ierr)
2062: c
2063: c  set up the preconditioner
2064: c
2065:       call PCDestroy(ctx(1),ierr)
2066:       call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
2067: CVAM  call PCSetType(ctx(1),PCJACOBI,ierr)
2068:       call PCSetType(ctx(1),PCLU,ierr)
2069: !      call PCSetVector(ctx(1),eold,ierr)
2070:       call PCSetOperators(ctx(1),ctx(2),ctx(2),
2071:      &              DIFFERENT_NONZERO_PATTERN,ierr)
2072:       call PCSetUp(ctx(1),ierr)

2074:       call VecDestroy(eold,ierr)


2077:       return
2078:       end
2079:       subroutine Setmat(ctx,ierr)

2081:       implicit none

2083:  #include finclude/petscsys.h
2084:  #include finclude/petscvec.h
2085:  #include finclude/petscmat.h
2086: !  Common blocks:
2087: #include "comd.h"
2088: #include "tube.h"

2090: !  Input/output variables:
2091:       PetscFortranAddr ctx(*)
2092:       integer  ierr

2094: !  Local variables:
2095:       PetscScalar      lx_v(0:1)
2096:       PetscOffset      lx_i

2098:       double precision xmult, himh, hiph, diag, upper, lower
2099:       double precision hi, hip1, him1
2100:       double precision
2101:      &             rhoee,  rhoe,  rhop,  rhow,  rhoww,
2102:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
2103:      &             ergee,  erge,  ergp,  ergw,  ergww,
2104:      &                     ue,    up,    uw
2105:       double precision see, sep, sew, seef, sewf, tef, twf,
2106:      &                 ref, rwf, kef, kwf, xmulte, xmultw
2107: c
2108:       integer  im, nx, ny
2109: c
2110: c     get pointers to xold
2111: c
2112:       call VecGetArray(ctx(4),lx_v,lx_i,ierr)
2113: 

2115: c
2116: c############################
2117: c
2118: c loop over all cells begin
2119: c
2120: c############################
2121: c
2122:       do im = 1,mx
2123: c
2124: c  set scalars
2125: c 
2126:          call Setpbc(im,lx_v(lx_i),
2127:      &             rhoee,  rhoe,  rhop,  rhow,  rhoww,
2128:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
2129:      &             ergee,  erge,  ergp,  ergw,  ergww,
2130:      &                     ue,    up,    uw)
2131: c
2132: c  set diffusion coefficients
2133: c
2134:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
2135:         sep = (ergp/rhop) - (0.5d+0 * up * up)
2136:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

2138:         seef = 0.5d+0 * (see + sep)
2139:         sewf = 0.5d+0 * (sew + sep)

2141:         tef = seef / csubv
2142:         twf = sewf / csubv

2144:         ref = 0.5d+0 * (rhoe + rhop)
2145:         rwf = 0.5d+0 * (rhow + rhop)

2147:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2148:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

2150:         if (wilson) then
2151:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2152:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2153:         endif
2154: c
2155: c  set coefficients
2156: c
2157:          xmult = dt / (dx * dx * csubv)

2159:          xmulte = xmult * kef
2160:          xmultw = xmult * kwf

2162:          upper = -(xmulte / rhoe)
2163:          lower = -(xmultw / rhow)

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

2167: c
2168: c  load coefficients into the matrix
2169: c
2170:          call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)

2172:          if (im .eq. 1) then
2173:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper,
2174:      1                       INSERT_VALUES,ierr)
2175:          elseif (im .eq. mx) then
2176:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2177:      1                       INSERT_VALUES,ierr)
2178:          else
2179:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper,
2180:      1                       INSERT_VALUES,ierr)
2181:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2182:      1                       INSERT_VALUES,ierr)
2183:          endif


2186:       enddo
2187: c
2188: c############################
2189: c
2190: c loop over all cells end
2191: c
2192: c############################
2193: c
2194: 
2195: c
2196: c  final load of matrix
2197: c
2198:       call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2199:       call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)

2201:       if (debug) then
2202:         call MatGetSize(ctx(2),nx,ny,ierr)
2203:         write(*,*) 'in setup nx = ',nx,' ny = ',ny
2204:         call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2205:       endif

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



2211:       return
2212:       end
2213:       subroutine Setpbc(i,x,
2214:      &             rhoee,  rhoe,  rhop,  rhow,  rhoww,
2215:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
2216:      &             ergee,  erge,  ergp,  ergw,  ergww,
2217:      &                     vele,  velp,  velw)

2219:       implicit none

2221: !  Common blocks:
2222: #include "comd.h"

2224: !  Input/output variables:
2225:       PetscScalar   x(mx*neq)
2226:       integer  i
2227:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2228:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2229:       double precision ergee,  erge,  ergp,  ergw,  ergww
2230:       double precision         vele,  velp,  velw

2232: !  Local variables:
2233:       integer  jr, jru, je

2235: c
2236: c  set pointers
2237: c
2238:       jr  = (neq*i) - 2
2239:       jru = (neq*i) - 1
2240:       je  = (neq*i)

2242:       if (debug) then
2243:         write(*,*)
2244:         write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2245:         write(*,*)
2246:       endif
2247: 
2248:       if (i .eq. 1) then

2250:         rhoee = x(jr+(2*neq))
2251:         rhoe  = x(jr+neq)
2252:         rhop  = x(jr)
2253:         rhow  = x(jr)
2254:         rhoww = x(jr)

2256:         rhouee = x(jru+(2*neq))
2257:         rhoue  = x(jru+neq)
2258:         rhoup  = x(jru)
2259:         rhouw  = x(jru)
2260:         rhouww = x(jru)

2262:         ergee = x(je+(2*neq))
2263:         erge  = x(je+neq)
2264:         ergp  = x(je)
2265:         ergw  = x(je)
2266:         ergww = x(je)

2268:         velw = 0.0d+0
2269:         velp = rhoup/rhop
2270:         vele = rhoue/rhoe

2272:       elseif (i .eq. 2) then

2274:         rhoee = x(jr+(2*neq))
2275:         rhoe  = x(jr+neq)
2276:         rhop  = x(jr)
2277:         rhow  = x(jr-neq)
2278:         rhoww = x(jr-neq)

2280:         rhouee = x(jru+(2*neq))
2281:         rhoue  = x(jru+neq)
2282:         rhoup  = x(jru)
2283:         rhouw  = x(jru-neq)
2284:         rhouww = x(jru-neq)

2286:         ergee = x(je+(2*neq))
2287:         erge  = x(je+neq)
2288:         ergp  = x(je)
2289:         ergw  = x(je-neq)
2290:         ergww = x(je-neq)

2292:         velw = rhouw/rhow
2293:         velp = rhoup/rhop
2294:         vele = rhoue/rhoe

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

2298:         rhoee = x(jr+neq)
2299:         rhoe  = x(jr+neq)
2300:         rhop  = x(jr)
2301:         rhow  = x(jr-neq)
2302:         rhoww = x(jr-(2*neq))

2304:         rhouee = x(jru+neq)
2305:         rhoue  = x(jru+neq)
2306:         rhoup  = x(jru)
2307:         rhouw  = x(jru-neq)
2308:         rhouww = x(jru-(2*neq))

2310:         ergee = x(je+neq)
2311:         erge  = x(je+neq)
2312:         ergp  = x(je)
2313:         ergw  = x(je-neq)
2314:         ergww = x(je-(2*neq))

2316:         velw = rhouw/rhow
2317:         velp = rhoup/rhop
2318:         vele = rhoue/rhoe

2320:       elseif (i .eq. mx) then

2322:         rhoee = x(jr)
2323:         rhoe  = x(jr)
2324:         rhop  = x(jr)
2325:         rhow  = x(jr-neq)
2326:         rhoww = x(jr-(2*neq))

2328:         rhouee = x(jru)
2329:         rhoue  = x(jru)
2330:         rhoup  = x(jru)
2331:         rhouw  = x(jru-neq)
2332:         rhouww = x(jru-(2*neq))

2334:         ergee = x(je)
2335:         erge  = x(je)
2336:         ergp  = x(je)
2337:         ergw  = x(je-neq)
2338:         ergww = x(je-(2*neq))

2340:         velw = rhouw/rhow
2341:         velp = rhoup/rhop
2342:         vele = 0.0d+0

2344:       else

2346:         rhoee = x(jr+(2*neq))
2347:         rhoe  = x(jr+neq)
2348:         rhop  = x(jr)
2349:         rhow  = x(jr-neq)
2350:         rhoww = x(jr-(2*neq))

2352:         rhouee = x(jru+(2*neq))
2353:         rhoue  = x(jru+neq)
2354:         rhoup  = x(jru)
2355:         rhouw  = x(jru-neq)
2356:         rhouww = x(jru-(2*neq))

2358:         ergee = x(je+(2*neq))
2359:         erge  = x(je+neq)
2360:         ergp  = x(je)
2361:         ergw  = x(je-neq)
2362:         ergww = x(je-(2*neq))

2364:         velw = rhouw/rhow
2365:         velp = rhoup/rhop
2366:         vele = rhoue/rhoe

2368:       endif

2370:       if (debug) then
2371:          write(*,*)
2372:          write(*,*) 'in Setpbc ',i,jr,jru,je
2373:          write(*,*) 'mx = ',mx
2374:          write(*,*)
2375:       endif


2378:       return
2379:       end
2380:       subroutine Setpbcn(rhoee,  rhoe,  rhop,  rhow,  rhoww,
2381:      &                   rhouee, rhoue, rhoup, rhouw, rhouww,
2382:      &                   ergee,  erge,  ergp,  ergw,  ergww,
2383:      &                           ue,    up,    uw,           jbc)

2385:       implicit none

2387: !  Common blocks:
2388: #include "comd.h"

2390: !  Input/output variables:
2391:       integer  jbc
2392:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2393:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2394:       double precision ergee,  erge,  ergp,  ergw,  ergww
2395:       double precision         ue,    up,    uw

2397: !  Local variables:

2399:       if (jbc .eq. 1) then
2400:          rhoww  = rhop
2401:          rhow   = rhop
2402:          rhouww = rhoup
2403:          rhouw  = rhoup
2404:          ergww  = ergp
2405:          ergw   = ergp
2406:          uw     = 0.0d+0
2407:       elseif  (jbc .eq. 2) then
2408:          rhoww  = rhow
2409:          rhouww = rhouw
2410:          ergww  = ergw
2411:          uw     = rhouw / rhow
2412:       else
2413:          uw = rhouw / rhow
2414:       endif

2416:       if (jbc .eq. mx) then
2417:          rhoee  = rhop
2418:          rhoe   = rhop
2419:          rhouee = rhoup
2420:          rhoue  = rhoup
2421:          ergee  = ergp
2422:          erge   = ergp
2423:          ue     = 0.0d+0
2424:       elseif (jbc .eq. mx-1) then
2425:          rhoee  = rhoe
2426:          rhouee = rhoue
2427:          ergee  = erge
2428:          ue     = rhoue / rhoe
2429:       else
2430:          ue     = rhoue / rhoe
2431:       endif

2433:       up = rhoup / rhop

2435:       if (debug) then
2436:          write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2437:       endif


2440:       return
2441:       end
2442:       double precision function cont
2443:      &            (rhoee,  rhoe,  rhop,  rhow,  rhoww,
2444:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
2445:      &             ergee,  erge,  ergp,  ergw,  ergww,
2446:      &                                      jcont,xold)
2447: c
2448: c  This function computes the residual
2449: c  for the 1-D continuity equation
2450: c
2451: c
2452:       implicit none

2454:       include 'comd.h'
2455:       include 'tube.h'
2456: c
2457: c     input variables
2458: c
2459:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2460:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2461:       double precision ergee,  erge,  ergp,  ergw,  ergww
2462:       double precision xold(mx*neq)
2463: c
2464:       integer jcont
2465: c
2466: c     local variables
2467: c
2468:       double precision theta1
2469:       integer jr
2470: c
2471: c  new
2472: c
2473:       double precision velfw, velfe
2474:       double precision vele,velp,velw
2475:       double precision fluxe, fluxw
2476:       double precision urhoe, urhow
2477:       double precision source
2478: c
2479: c old
2480: c
2481:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2482:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2483:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2484:       double precision teoee, teoe, teop, teow, teoww,
2485:      &                 uoe, uoee, uop, uow, uoww
2486:       double precision velfow, velfoe
2487:       double precision veloe,velop,velow
2488:       double precision fluxoe, fluxow
2489:       double precision urhooe, urhoow
2490:       double precision sourceo
2491: c
2492: c functions
2493: c
2494:       double precision godunov2
2495:       double precision upwind, fluxlim
2496: c
2497: c
2498: c ******************************************************************
2499: c
2500: c
2501: c
2502:       if (debug) then
2503:         write(*,*)
2504:         write(*,*) 'in cont',jcont,' ihod = ',ihod
2505:         write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2506:         write(*,*) 'rhop = ',rhop
2507:         write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2508:         write(*,*)
2509:       endif

2511:       jr = (neq*jcont) - 2

2513: c########################
2514: c
2515: c      NEW
2516: c
2517: c########################

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

2524:       velfe = 0.5d+0 * (vele + velp)
2525:       velfw = 0.5d+0 * (velw + velp)

2527:       if (ihod .eq. 1) then

2529:         urhoe = upwind(rhop,rhoe,velfe)
2530:         urhow = upwind(rhow,rhop,velfw)

2532:       elseif (ihod .eq. 2) then

2534:         urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2535:         urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)

2537:       endif

2539:       if (ihod .eq. 3) then
2540:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2541:      &                             rhouw,rhoup,rhoue,rhouee,
2542:      &                             ergw, ergp, erge, ergee,1)
2543:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2544:      &                             rhouww,rhouw,rhoup,rhoue,
2545:      &                             ergww, ergw, ergp, erge,1)
2546:       else
2547:         fluxe = (dt/dx) * urhoe
2548:         fluxw = (dt/dx) * urhow
2549:       endif


2552:       source = 0.0d+0

2554: c########################
2555: c
2556: c      OLD
2557: c
2558: c########################

2560:       call Setpbc(jcont,xold,
2561:      &             rhooee,  rhooe,  rhoop,  rhoow,  rhooww,
2562:      &             rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2563:      &             ergoee,  ergoe,  ergop,  ergow,  ergoww,
2564:      &                      veloe,  velop,  velow)

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

2571:       velfoe = 0.5d+0 * (veloe + velop)
2572:       velfow = 0.5d+0 * (velow + velop)


2575:       if (ihod .eq. 1) then

2577:         urhooe = upwind(rhoop,rhooe,velfoe)
2578:         urhoow = upwind(rhoow,rhoop,velfow)

2580:       elseif (ihod .eq. 2) then

2582:         urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2583:         urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)

2585:       endif

2587:       if (ihod .eq. 3) then
2588:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
2589:      &                              rhouow,rhouop,rhouoe,rhouoee,
2590:      &                              ergow, ergop, ergoe, ergoee,1)
2591:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
2592:      &                              rhouoww,rhouow,rhouop,rhouoe,
2593:      &                              ergoww, ergow, ergop, ergoe,1)
2594:       else
2595:         fluxoe = (dt/dx) * urhooe
2596:         fluxow = (dt/dx) * urhoow
2597:       endif

2599:       sourceo = 0.0d+0


2602: c########################
2603: c
2604: c      FUNCTION
2605: c
2606: c########################

2608:       theta1 = 1.0d+0 - theta
2609:       cont =  (rhop - xold(jr))
2610:      &    + (  theta  * ( (fluxe  - fluxw ) - source  )  )
2611:      &    + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
2612: CVAM
2613:       if (probnum .eq. 3) then
2614:         cont = 0.0d+0
2615:       endif
2616: CVAM


2619:       if (debug) then
2620:        write(*,*)
2621:        write(*,*) 'cont(',jcont,') = ',cont
2622:        write(*,*) 'theta = ',theta,'rhop = ',rhop
2623:        write(*,*) 'source = ',source,' sourceo = ',sourceo
2624:        write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2625:        write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2626:        write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2627:        write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2628:        write(*,*)
2629:       endif

2631:       return
2632:       end
2633:       double precision function  eexact(x,t)

2635:       implicit none

2637:       double precision x,t
2638:       double precision xot, head, tail, contact, ufan
2639:       double precision xpow, grat, urat
2640:       double precision uexact


2643:       logical debug

2645:       include 'tube.h'

2647:       debug = .false.


2650:       if (t .le. 0.0d+0) then
2651:         if (x .gt. 0.0d+0) then
2652:           eexact = e1
2653:         else
2654:           eexact = e4
2655:         endif
2656:       else

2658:        xot = x/t
2659:        head = -a4
2660:        tail = v3 - a3
2661:        contact = v2

2663:        if (xot .lt. head) then
2664:           eexact = e4
2665:        elseif (xot .gt. sspd) then
2666:           eexact = e1
2667:        elseif (xot .gt. contact) then
2668:           eexact = e2
2669:        elseif (xot .gt. tail) then
2670:           eexact = e3
2671:        else
2672:           ufan = uexact(x,t)
2673:           grat = (gamma - 1.0d+0) / 2.0d+0
2674:           xpow = 2.0d+0
2675:           urat = ufan / a4
2676:           eexact = e4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
2677:        endif

2679:       endif


2682:       if (debug) then
2683:         write(*,*)
2684:         write(*,*) 'eexact(',x,',',t,') = ',eexact
2685:         write(*,*)
2686:       endif

2688:       return
2689:       end
2690:       subroutine eigen(ht,uht)
2691: c23456789012345678901234567890123456789012345678901234567890123456789012
2692: c
2693: c          subroutine eigen
2694: c
2695: c  This subroutine computes the eigen values and eigen vectors
2696: c
2697: c23456789012345678901234567890123456789012345678901234567890123456789012


2700: c#######################################################################

2702:       implicit none

2704:       include 'comd.h'

2706:       double precision ht, uht

2708:       double precision ut, at, lam1, lam2


2711: c#######################################################################

2713:       ut = uht / ht
2714:       at = sqrt( ht)

2716:       lam1 = ut - at
2717:       lam2 = ut + at

2719:       eigval(1) = lam1
2720:       eigval(2) = lam2

2722:       eigvec(1,1) = 1.0d+0
2723:       eigvec(2,1) = lam1
2724:       eigvec(1,2) = 1.0d+0
2725:       eigvec(2,2) = lam2

2727:       rinv(1,1) =  lam2 / (2.0d+0 * at)
2728:       rinv(2,1) = -lam1 / (2.0d+0 * at)
2729:       rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2730:       rinv(2,2) =  1.0d+0 / (2.0d+0 * at)


2733:       return
2734:       end
2735:       subroutine eigene(r,ru,e,l1,l2,l3)
2736: c23456789012345678901234567890123456789012345678901234567890123456789012
2737: c
2738: c          subroutine eigene
2739: c
2740: c  This subroutine computes the eigen values for the entropy fix
2741: c
2742: c23456789012345678901234567890123456789012345678901234567890123456789012


2745: c#######################################################################

2747:       implicit none

2749:       include 'tube.h'

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

2753:       double precision p,u,a

2755:       double precision eos
2756:       integer ierr

2758:       logical debug


2761: c#######################################################################

2763:       debug = .false.

2765:       if (debug) then
2766:          write(*,*)
2767:          write(*,*) 'gamma = ',gamma
2768:          write(*,*) 'r,ru,e = ',r,ru,e
2769:          write(*,*)
2770:       endif

2772:       p = eos(r,ru,e)
2773:       u = ru/r
2774:       if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2775:          write(*,*)
2776:          write(*,*) 'gamma = ',gamma
2777:          write(*,*) 'r = ',r
2778:          write(*,*) 'p = ',p
2779:          write(*,*)
2780:          call PetscFinalize(ierr)
2781:          stop
2782:       endif
2783:       a = sqrt((gamma * p)/r)

2785:       if (debug) then
2786:          write(*,*)
2787:          write(*,*) 'p,u,a = ',p,u,a
2788:          write(*,*)
2789:       endif

2791:       l1 = u - a
2792:       l2 = u
2793:       l3 = u + a

2795:       return
2796:       end
2797:       double precision function energy
2798:      &               (rhoee,  rhoe,  rhop,  rhow,  rhoww,
2799:      &                rhouee, rhoue, rhoup, rhouw, rhouww,
2800:      &                ergee,  erge,  ergp,  ergw,  ergww,
2801:      &                                      jerg,xold)
2802: c
2803: c  This function computes the residual
2804: c  for the 1-D energy equation
2805: c
2806: c
2807:       implicit none

2809:       include 'comd.h'
2810:       include 'tube.h'
2811: c
2812: c     input variables
2813: c
2814:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2815:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2816:       double precision ergee,  erge,  ergp,  ergw,  ergww
2817:       double precision xold(mx*neq)
2818: c
2819:       integer jerg
2820: c
2821: c     local variables
2822: c
2823:       double precision theta1
2824:       integer je
2825: c
2826: c  new
2827: c
2828:       double precision velfw, velfe
2829:       double precision vele,velp,velw
2830:       double precision fluxe, fluxw
2831:       double precision uepe, uepw
2832:       double precision ue, up, uw
2833:       double precision see, sep, sew
2834:       double precision seef, sewf
2835:       double precision upe, upw
2836:       double precision presse, pressw
2837:       double precision source
2838:       double precision te, tp, tw
2839:       double precision tef, twf, ref, rwf
2840:       double precision kef, kwf
2841:       double precision hflxe, hflxw
2842: c
2843: c old
2844: c
2845:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2846:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2847:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2848:       double precision velfow, velfoe
2849:       double precision veloe,velop,velow
2850:       double precision fluxoe, fluxow
2851:       double precision uepoe, uepow
2852:       double precision uoe, uop, uow
2853:       double precision seoe, seop, seow
2854:       double precision seoef, seowf
2855:       double precision upoe, upow
2856:       double precision pressoe, pressow
2857:       double precision sourceo
2858:       double precision toe, top, tow
2859:       double precision toef, towf, roef, rowf
2860:       double precision koef, kowf
2861:       double precision hflxoe, hflxow
2862: c
2863: c functions
2864: c
2865:       double precision godunov2, eos
2866:       double precision upwind, fluxlim

2868: c
2869: c
2870: c ******************************************************************
2871: c
2872: c
2873: c
2874:       je = (neq*jerg)

2876: c########################
2877: c
2878: c      NEW
2879: c
2880: c########################

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

2887:       pressw  = eos(rhow, rhouw, ergw)
2888:       presse  = eos(rhoe, rhoue, erge)

2890:       uw = rhouw / rhow
2891:       up = rhoup / rhop
2892:       ue = rhoue / rhoe

2894:       upw = uw * pressw
2895:       upe = ue * presse

2897:       velfe = 0.5d+0 * (vele + velp)
2898:       velfw = 0.5d+0 * (velw + velp)

2900:       if (ihod .eq. 1) then

2902:         uepe = upwind(ergp,erge,velfe)
2903:         uepw = upwind(ergw,ergp,velfw)

2905:       elseif (ihod .eq. 2) then

2907:         uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2908:         uepw = fluxlim(ergww,ergw,ergp,erge,velfw)

2910:       endif

2912:       if (ihod .eq. 3) then
2913:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2914:      &                             rhouw,rhoup,rhoue,rhouee,
2915:      &                             ergw, ergp, erge, ergee,3)
2916:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2917:      &                             rhouww,rhouw,rhoup,rhoue,
2918:      &                             ergww, ergw, ergp, erge,3)
2919:       else
2920:         fluxe = (dt/dx) * ( uepe  + (0.5d+0*upe) )
2921:         fluxw = (dt/dx) * ( uepw  + (0.5d+0*upw) )
2922:       endif
2923: c
2924: c  radiation
2925: c
2926:       if (kappa0 .eq. 0.0d+0) then
2927:         source = 0.0d+0
2928:       else

2930:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
2931:         sep = (ergp/rhop) - (0.5d+0 * up * up)
2932:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

2934:         seef = 0.5d+0 * (see + sep)
2935:         sewf = 0.5d+0 * (sew + sep)

2937:         te  = see / csubv
2938:         tp  = sep / csubv
2939:         tw  = sew / csubv

2941:         tef = seef / csubv
2942:         twf = sewf / csubv

2944:         ref = 0.5d+0 * (rhoe + rhop)
2945:         rwf = 0.5d+0 * (rhow + rhop)

2947:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2948:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

2950:         if (wilson) then
2951:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2952:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2953:         endif

2955:         if ( debug .and. (kef .gt. 1.0d+10) ) then
2956:           write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2957:         endif
2958:         if ( debug .and. (kwf .gt. 1.0d+10) ) then
2959:           write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2960:         endif

2962:         hflxe = kef * (te - tp) / dx
2963:         hflxw = kwf * (tp - tw) / dx

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

2967:       endif

2969: c########################
2970: c
2971: c      OLD
2972: c
2973: c########################

2975:       call Setpbc(jerg,xold,
2976:      &             rhooee,  rhooe,  rhoop,  rhoow,  rhooww,
2977:      &             rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2978:      &             ergoee,  ergoe,  ergop,  ergow,  ergoww,
2979:      &                      veloe,  velop,  velow)

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

2986:       pressow  = eos(rhoow, rhouow, ergow)
2987:       pressoe  = eos(rhooe, rhouoe, ergoe)

2989:       uow = rhouow / rhoow
2990:       uop = rhouop / rhoop
2991:       uoe = rhouoe / rhooe

2993:       upow = uow * pressow
2994:       upoe = uoe * pressoe

2996:       velfoe = 0.5d+0 * (veloe + velop)
2997:       velfow = 0.5d+0 * (velow + velop)


3000:       if (ihod .eq. 1) then

3002:         uepoe = upwind(ergop,ergoe,velfoe)
3003:         uepow = upwind(ergow,ergop,velfow)

3005:       elseif (ihod .eq. 2) then

3007:         uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
3008:         uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)

3010:       endif

3012:       if (ihod .eq. 3) then
3013:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3014:      &                              rhouow,rhouop,rhouoe,rhouoee,
3015:      &                              ergow, ergop, ergoe, ergoee,3)
3016:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3017:      &                              rhouoww,rhouow,rhouop,rhouoe,
3018:      &                              ergoww, ergow, ergop, ergoe,3)
3019:       else
3020:         fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
3021:         fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
3022:       endif

3024: c
3025: c  old radiation
3026: c
3027:       if (kappa0 .eq. 0.0d+0) then
3028:         sourceo = 0.0d+0
3029:       else

3031:         seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
3032:         seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
3033:         seow = (ergow/rhoow) - (0.5d+0 * uow * uow)

3035:         seoef = 0.5d+0 * (seoe + seop)
3036:         seowf = 0.5d+0 * (seow + seop)

3038:         toe  = seoe / csubv
3039:         top  = seop / csubv
3040:         tow  = seow / csubv

3042:         toef = seoef / csubv
3043:         towf = seowf / csubv

3045:         roef = 0.5d+0 * (rhooe + rhoop)
3046:         rowf = 0.5d+0 * (rhoow + rhoop)

3048:         koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
3049:         kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)

3051:         if (wilson) then
3052:            koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
3053:            kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
3054:         endif

3056:         if ( debug .and. (koef .gt. 1.0d+10) ) then
3057:           write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
3058:         endif
3059:         if ( debug .and. (kowf .gt. 1.0d+10) ) then
3060:           write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
3061:         endif

3063:         hflxoe = koef * (toe - top) / dx
3064:         hflxow = kowf * (top - tow) / dx

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

3068:       endif


3071: c########################
3072: c
3073: c      FUNCTION
3074: c
3075: c########################

3077: CVAM
3078:       if (probnum .eq. 3) then
3079:         fluxe  = 0.0d+0
3080:         fluxw  = 0.0d+0
3081:         fluxoe = 0.0d+0
3082:         fluxow = 0.0d+0
3083:       endif
3084: CVAM

3086:       theta1 = 1.0d+0 - theta
3087:       energy =  (ergp - xold(je))
3088:      &    + (  theta  * ( (fluxe  - fluxw ) - source  )  )
3089:      &    + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )

3091:       if (debug) then
3092:         write(*,*)
3093:         write(*,*) 'energy(',jerg,') = ',energy
3094:         write(*,*)
3095:         write(*,*) fluxe,fluxw
3096:         write(*,*) fluxoe,fluxow
3097:         write(*,*) source,sourceo
3098:         write(*,*)
3099:       endif

3101:       return
3102:       end
3103:       double precision function eos(r,ru,e)

3105:       implicit none

3107:       double precision r,ru,e

3109:       double precision se, u

3111:       integer ierr

3113:       logical debug

3115:       include "tube.h"

3117:       debug = .false.

3119:       if (debug) then
3120:         write(*,*)
3121:         write(*,*) 'in eos r,ru,e'
3122:         write(*,*) r,ru,e
3123:         write(*,*)
3124:       endif

3126:       u = ru/r

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

3131:       if (eos .lt. 0.0d+0) then
3132:          write(*,*)
3133:          write(*,*) 'eos = ',eos
3134:          write(*,*) 'gamma = ',gamma
3135:          write(*,*) 'r = ',r
3136:          write(*,*) 'se = ',se
3137:          write(*,*) 'e = ',e
3138:          write(*,*) 'u = ',u
3139:          write(*,*) 'ru = ',ru
3140:          call PetscFinalize(ierr)
3141:          write(*,*)
3142:          stop
3143:       endif

3145:       if (debug) then
3146:         write(*,*)
3147:         write(*,*) 'in eos u,se,eos'
3148:         write(*,*) u,se,eos
3149:         write(*,*)
3150:       endif


3153:       return
3154:       end
3155:       subroutine eval2

3157:       implicit none

3159:       double precision prat, grat, xnum, xdenom


3162:       logical debug

3164:       include 'tube.h'

3166:       debug = .false.

3168:       prat = p2/p1
3169:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

3171:       xnum = grat + prat
3172:       xdenom = 1.0d+0 + (prat * grat)
3173: 
3174:       e2 = e1 * prat * (xnum/xdenom)
3175: 


3178:       if (debug) then
3179:         write(*,*)
3180:         write(*,*) 'e1  = ',e1
3181:         write(*,*) 'e2  = ',e2
3182:       endif

3184:       return
3185:       end
3186:       subroutine exact0

3188:       implicit none

3190:       double precision tol, xn
3191:       double precision shockp, fprime

3193:       integer maxnewt, niter

3195:       logical found, debug

3197:       include 'tube.h'

3199:       debug = .false.

3201:       tol = 1.0d-10

3203:       maxnewt = 40
3204: 
3205:       a1 = sqrt(gamma*p1/r1)
3206:       a4 = sqrt(gamma*p4/r4)



3210:       found = .false.
3211:       niter = 0

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

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

3217:         niter = niter + 1

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

3221:         if (debug) then
3222:           write(*,*) niter,shockp(xn),xn
3223:         endif

3225:         if ( abs(shockp(xn)) .lt. tol ) then
3226:            found = .true.
3227:         endif

3229:         goto 10

3231:       endif

3233:       if (.not. found) then

3235:          write(*,*) 'newton failed'
3236:          write(*,*) xn,shockp(xn)
3237:          stop

3239:       endif

3241:       p2 = xn


3244:       if (debug) then
3245:         write(*,*)
3246:         write(*,*) 'p1  = ',p1
3247:         write(*,*) 'p2  = ',p2
3248:         write(*,*) 'p4  = ',p4
3249:         write(*,*)
3250:       endif

3252:       return
3253:       end
3254:       double precision function flux(r,ru,e,eqn)
3255: c23456789012345678901234567890123456789012345678901234567890123456789012
3256: c
3257: c          function flux
3258: c
3259: c  This function computes the flux at a face
3260: c
3261: c23456789012345678901234567890123456789012345678901234567890123456789012


3264: c#######################################################################

3266:       implicit none

3268:       include 'comd.h'
3269:       include 'tube.h'

3271:       double precision r, ru, e

3273:       integer eqn

3275:       double precision p,u

3277:       double precision eos


3280: c#######################################################################

3282:       p = eos(r,ru,e)
3283:       u = ru/r

3285:       if (eqn .eq. 1) then
3286:          flux = ru
3287:       elseif (eqn .eq. 2) then
3288:          flux = (u * ru) + p
3289:       else
3290:          flux = u * (e + p)
3291:       endif

3293:       return
3294:       end
3295:       double precision function fluxlim(fww,fw,fe,fee,vp)
3296: c23456789012345678901234567890123456789012345678901234567890123456789012
3297: c
3298: c          function fluxlim
3299: c
3300: c  this function computes the flux limited quick face value
3301: c
3302: c23456789012345678901234567890123456789012345678901234567890123456789012


3305: c#######################################################################

3307:       implicit none

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

3311:       double precision fd, fc, fu

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

3315:       double precision med, quick
3316: 
3317:       logical limit

3319: c#######################################################################

3321:       limit = .true.

3323:       if (vp .gt. 0.0d+0) then
3324:         fd = fe
3325:         fc = fw
3326:         fu = fww
3327:       else
3328:         fd = fw
3329:         fc = fe
3330:         fu = fee
3331:       endif

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

3335:       if (limit) then

3337:         beta = 0.25d+0
3338:         flc = 4.0d+0

3340:         f1 = fc
3341:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3342:         f3 = fu + ( flc * (fc - fu) )
3343:         f4 = med(f1,f2,f3)
3344:         fluxlim = vp * med(f1,f4,fhod)

3346:       else

3348:         fluxlim = vp * fhod

3350:       endif

3352:       return
3353:       end
3354:       double precision function fluxlim2(fww,fw,fe,fee,vp)
3355: c23456789012345678901234567890123456789012345678901234567890123456789012
3356: c
3357: c          function fluxlim2
3358: c
3359: c  this function computes the flux limited quick face value
3360: c
3361: c23456789012345678901234567890123456789012345678901234567890123456789012


3364: c#######################################################################

3366:       implicit none

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

3370:       double precision fd, fc, fu

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

3374:       double precision med, quick
3375: 
3376:       logical limit, debug

3378: c#######################################################################

3380:       debug = .false.

3382:       if (debug) then
3383:         write(*,*)
3384:         write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3385:         write(*,*) fee,fe,fw,fww
3386:         write(*,*)
3387:       endif

3389:       limit = .true.

3391:       if (vp .gt. 0.0d+0) then
3392:         fd = fe
3393:         fc = fw
3394:         fu = fww
3395:       else
3396:         fd = fw
3397:         fc = fe
3398:         fu = fee
3399:       endif

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

3403:       if (limit) then

3405:         beta = 0.25d+0
3406:         flc = 4.0d+0

3408:         f1 = fc
3409:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3410:         f3 = fu + ( flc * (fc - fu) )
3411:         f4 = med(f1,f2,f3)
3412:         fluxlim2 =  med(f1,f4,fhod)

3414:       else

3416:         fluxlim2 = fhod

3418:       endif

3420:       return
3421:       end
3422:       double precision function fprime(x)

3424:       implicit none

3426:       double precision  x, eps
3427:       double precision  shockp

3429:       eps = 1.0d-8

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

3433:       return
3434:       end
3435:       double precision function godent(uhl, uhr, hl, hr, eqn)
3436: c23456789012345678901234567890123456789012345678901234567890123456789012
3437: c
3438: c          function godent
3439: c
3440: c  this function computes the roe/godunov face value plus entropy fix
3441: c
3442: c23456789012345678901234567890123456789012345678901234567890123456789012


3445: c#######################################################################

3447:       implicit none

3449:       include 'comd.h'

3451:       double precision uhl, uhr, hl, hr, ht, uht
3452:       double precision lamr1, lamr2, laml1, laml2
3453:       double precision deltal1, deltal2

3455:       integer eqn

3457:       double precision sum

3459:       double precision flux

3461:       integer i, j

3463: 

3465: c#######################################################################

3467:       if (debug) then
3468:         write(*,*) 'in godent eqn = ',eqn
3469:       endif

3471: c      do i = 1,neq
3472: c        fr(i) = flux(uhr,hr,i)
3473: c        fl(i) = flux(uhl,hl,i)
3474: c      enddo

3476:       deltau(1) = hr - hl
3477:       deltau(2) = uhr - uhl

3479:       call roestat(uhl, uhr, hl,hr,ht, uht)

3481:       call eigen(ht,uht)

3483:       do i = 1,neq
3484:         sum = 0.0d+0
3485:         do j = 1,neq
3486:           sum = sum + ( rinv(i,j) * deltau(j) )
3487:         enddo
3488:         alpha(i) = sum
3489:       enddo

3491:       deltal1 = 0.0d+0
3492:       deltal2 = 0.0d+0

3494: c      call eigene(hr,uhr,lamr1, lamr2)
3495: c      call eigene(hl,uhl,laml1, laml2)
3496: c
3497: c  1st eigen
3498: c
3499:       if ( (laml1 .lt. 0.0d+0) .and.
3500:      &     (lamr1 .gt. 0.0d+0) ) then

3502: CVAM     deltal1 = 4.0d+0 * (lamr1 - laml1)
3503:          deltal1 = 4.0d+0 * (lamr1 - laml1) + 1.0d-2

3505:          if ( abs(eigval(1)) .lt. (0.5d+0 * deltal1) ) then
3506:              eigval(1) = (  ( (eigval(1) ** 2) / deltal1 )
3507:      &                    + ( 0.25d+0 * deltal1 )  )
3508:          endif

3510:       endif
3511: c
3512: c  2nd eigen
3513: c
3514:       if ( (laml2 .lt. 0.0d+0) .and.
3515:      &     (lamr2 .gt. 0.0d+0) ) then

3517:          deltal2 = 4.0d+0 * (lamr2 - laml2)

3519:          if ( abs(eigval(2)) .lt. (0.5d+0 * deltal2) ) then
3520:              eigval(2) = (  ( (eigval(2) ** 2) / deltal2 )
3521:      &                    + ( 0.25d+0 * deltal2 )  )
3522:          endif

3524:       endif

3526:       if (debug) then
3527:          write(*,*)
3528:          write(*,*) 'godent debug'
3529:          write(*,*) laml1, laml2, lamr1, lamr2
3530:          write(*,*) deltal1, deltal2
3531:          write(*,*)
3532:       endif

3534:       do i = 1,neq
3535:         sum = 0.0d+0
3536:         do j = 1,neq
3537:           sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3538:         enddo
3539:         xnumdif(i) = sum
3540:       enddo

3542:       do i = 1,neq
3543:         froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3544:       enddo

3546:       godent = froe(eqn)

3548:       return
3549:       end
3550:       double precision function godunov(uhl, uhr, hl, hr, eqn)
3551: c23456789012345678901234567890123456789012345678901234567890123456789012
3552: c
3553: c          function godunov
3554: c
3555: c  this function computes the roe/godunov face value
3556: c
3557: c23456789012345678901234567890123456789012345678901234567890123456789012


3560: c#######################################################################

3562:       implicit none

3564:       include 'comd.h'

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

3568:       integer eqn

3570:       double precision sum

3572:       double precision flux

3574:       integer i, j
3575: 

3577: c#######################################################################

3579: c      do i = 1,neq
3580: c        fr(i) = flux(uhr,hr,i)
3581: c        fl(i) = flux(uhl,hl,i)
3582: c      enddo

3584:       deltau(1) = hr - hl
3585:       deltau(2) = uhr - uhl

3587:       call roestat(uhl, uhr, hl,hr,ht, uht)

3589:       call eigen(ht,uht)

3591:       do i = 1,neq
3592:         sum = 0.0d+0
3593:         do j = 1,neq
3594:           sum = sum + ( rinv(i,j) * deltau(j) )
3595:         enddo
3596:         alpha(i) = sum
3597:       enddo

3599:       do i = 1,neq
3600:         sum = 0.0d+0
3601:         do j = 1,neq
3602:           sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3603:         enddo
3604:         xnumdif(i) = sum
3605:       enddo

3607:       do i = 1,neq
3608:         froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3609:       enddo

3611:       godunov = froe(eqn)


3614:       return
3615:       end
3616:       double precision function godunov2
3617:      &            (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3618: c23456789012345678901234567890123456789012345678901234567890123456789012
3619: c
3620: c          function godunov2
3621: c
3622: c  this function computes the roe/godunov2 face value
3623: c
3624: c23456789012345678901234567890123456789012345678901234567890123456789012


3627: c#######################################################################

3629:       implicit none

3631:       include 'comd.h'
3632:       include 'tube.h'

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

3636:       integer eqn

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

3640:       double precision godunov, godent, hlle


3643: c#######################################################################

3645:       if (gorder .eq. 1) then
3646:         rrg  = rr
3647:         rlg  = rl
3648:         rurg = rur
3649:         rulg = rul
3650:         erg  = er
3651:         elg  = el
3652:       else
3653:         call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
3654:      &               rrg, rlg,rurg, rulg, erg, elg)
3655:       endif

3657: CVAM  if (ientro .eq. 0) then
3658: CVAM     godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3659: CVAM  elseif(ientro .eq. 1) then
3660: CVAM     godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3661: CVAM  else
3662:          godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3663: CVAM  endif


3666:       return
3667:       end
3668:       double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3669: c23456789012345678901234567890123456789012345678901234567890123456789012
3670: c
3671: c          function hlle
3672: c
3673: c  this function computes the roe/hlle face value
3674: c
3675: c23456789012345678901234567890123456789012345678901234567890123456789012


3678: c#######################################################################

3680:       implicit none

3682:       include 'comd.h'
3683:       include 'tube.h'

3685:       double precision rrg,rlg,rurg,rulg,erg,elg
3686:       integer eqn

3688:       double precision laml1, laml2, laml3
3689:       double precision lamr1, lamr2, lamr3
3690:       double precision sl, sr


3693:       double precision flux

3695:       integer i, j, ispeed
3696: 

3698: c#######################################################################

3700:       ispeed = 1

3702:       do i = 1,neq
3703:         fr(i) = flux(rrg,rurg,erg,i)
3704:         fl(i) = flux(rlg,rulg,elg,i)
3705:       enddo

3707:       deltau(1) = rrg  - rlg
3708:       deltau(2) = rurg - rulg
3709:       deltau(3) = erg  - elg

3711: CVAM  call roestat(uhl,uhr,hl,hr,ht,uht)

3713: CVAM  call eigene(ht,uht,lamt1, lamt2)
3714:       call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3715:       call eigene(rlg,rulg,elg,laml1,laml2,laml3)

3717: CVAM  if (ispeed .eq. 1) then
3718: CVAM    sl = min(laml1,lamt1)
3719: CVAM    sr = max(lamt2,lamr2)
3720: CVAM  else
3721:         sl = min(laml1,lamr1)
3722:         sr = max(laml3,lamr3)
3723: CVAM  endif


3726:       do i = 1,neq
3727:         froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )
3728:      &            / (sr-sl)
3729:       enddo

3731:       hlle = froe(eqn)


3734:       return
3735:       end
3736:       double precision function med(x1,x2,x3)
3737: c23456789012345678901234567890123456789012345678901234567890123456789012
3738: c
3739: c          function med
3740: c
3741: c  this function computes the median of three numbers
3742: c
3743: c23456789012345678901234567890123456789012345678901234567890123456789012


3746: c#######################################################################

3748:       implicit none

3750:       double precision x1, x2, x3
3751:       double precision xhi, xlo

3753: c#######################################################################

3755:       xhi = max(x1,x2,x3)
3756:       xlo = min(x1,x2,x3)

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

3760:       return
3761:       end
3762:       double precision function mom
3763:      &            (rhoee,  rhoe,  rhop,  rhow,  rhoww,
3764:      &             rhouee, rhoue, rhoup, rhouw, rhouww,
3765:      &             ergee,  erge,  ergp,  ergw,  ergww,
3766:      &                                      jmom,xold)
3767: c
3768: c  This function computes the residual
3769: c  for the 1-D momentum equation
3770: c
3771: c
3772:       implicit none

3774:       include 'comd.h'
3775:       include 'tube.h'
3776: c
3777: c     input variables
3778: c
3779:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
3780:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
3781:       double precision ergee,  erge,  ergp,  ergw,  ergww
3782:       double precision xold(mx*neq)
3783: c
3784:       integer jmom
3785: c
3786: c     local variables
3787: c
3788:       double precision theta1
3789:       integer jru
3790: c
3791: c  new
3792: c
3793:       double precision velfw, velfe
3794:       double precision vele,velp,velw
3795:       double precision fluxe, fluxw
3796:       double precision uurhoe, uurhow
3797:       double precision pressee, presse, pressp,pressw, pressww
3798:       double precision rupee, rupe, rupp, rupw, rupww
3799:       double precision uee, ue, up, uw, uww
3800:       double precision source
3801: c
3802: c old
3803: c
3804:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
3805:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3806:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
3807:       double precision velfow, velfoe
3808:       double precision veloe,velop,velow
3809:       double precision fluxoe, fluxow
3810:       double precision uurhooe, uurhoow
3811:       double precision pressoee, pressoe, pressop, pressow, pressoww
3812:       double precision rupoee, rupoe, rupop, rupow, rupoww
3813:       double precision uoee, uoe, uop, uow, uoww
3814:       double precision sourceo

3816:       double precision eps
3817: c
3818: c functions
3819: c
3820:       double precision godunov2, eos
3821:       double precision upwind, fluxlim
3822: c
3823: c
3824: c ******************************************************************
3825: c
3826: c
3827:       eps = 1.0d-32
3828: c
3829:       jru = (neq*jmom) - 1

3831: c########################
3832: c
3833: c      NEW
3834: c
3835: c########################

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

3842:       presse  = eos(rhoe, rhoue, erge )
3843:       pressw  = eos(rhow, rhouw, ergw )

3845:       velfe = 0.5d+0 * (vele + velp)
3846:       velfw = 0.5d+0 * (velw + velp)

3848:       if (ihod .eq. 1) then

3850:         uurhoe = upwind(rhoup,rhoue,velfe)
3851:         uurhow = upwind(rhouw,rhoup,velfw)

3853:       elseif (ihod .eq. 2) then

3855:         uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3856:         uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)

3858:       endif

3860:       if (ihod .eq. 3) then
3861:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
3862:      &                             rhouw,rhoup,rhoue,rhouee,
3863:      &                             ergw, ergp, erge, ergee,2)
3864:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
3865:      &                             rhouww,rhouw,rhoup,rhoue,
3866:      &                             ergww, ergw, ergp, erge,2)
3867:       else
3868:         fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3869:         fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3870:       endif


3873:       source = 0.0d+0

3875: c########################
3876: c
3877: c      OLD
3878: c
3879: c########################

3881:       call Setpbc(jmom,xold,
3882:      &             rhooee,  rhooe,  rhoop,  rhoow,  rhooww,
3883:      &             rhouoee, rhouoe, rhouop, rhouow, rhouoww,
3884:      &             ergoee,  ergoe,  ergop,  ergow,  ergoww,
3885:      &                      veloe,  velop,  velow)

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

3892:       pressoe  = eos(rhooe, rhouoe, ergoe)
3893:       pressow  = eos(rhoow, rhouow, ergow)

3895:       velfoe = 0.5d+0 * (veloe + velop)
3896:       velfow = 0.5d+0 * (velow + velop)

3898:       if (ihod .eq. 1) then

3900:         uurhooe = upwind(rhouop,rhouoe,velfoe)
3901:         uurhoow = upwind(rhouow,rhouop,velfow)

3903:       elseif (ihod .eq. 2) then

3905:         uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3906:         uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)

3908:       endif

3910:       if (ihod .eq. 3) then
3911:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3912:      &                              rhouow,rhouop,rhouoe,rhouoee,
3913:      &                              ergow, ergop, ergoe, ergoee,2)
3914:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3915:      &                              rhouoww,rhouow,rhouop,rhouoe,
3916:      &                              ergoww, ergow, ergop, ergoe,2)
3917:       else
3918:         fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3919:         fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3920:       endif

3922:       sourceo = 0.0d+0


3925: c########################
3926: c
3927: c      FUNCTION
3928: c
3929: c########################

3931:       theta1 = 1.0d+0 - theta
3932:       mom =  (rhoup - xold(jru))
3933:      &    + (  theta  * ( (fluxe  - fluxw ) - source  )  )
3934:      &    + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
3935: CVAM
3936:       if (probnum .eq. 3) then
3937:         mom = 0.0d+0
3938:       endif
3939: CVAM
3940:       if (debug) then
3941:         write(*,*)
3942:         write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3943:         write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3944:         write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3945:         write(*,*) 'presse = ',presse,'pressw = ',pressw
3946:         write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3947:         write(*,*)
3948:       endif

3950:       return
3951:       end
3952:       double precision function quick(fd, fc, fu)
3953: c23456789012345678901234567890123456789012345678901234567890123456789012
3954: c
3955: c          function quick
3956: c
3957: c  this function computes the quick face value
3958: c
3959: c23456789012345678901234567890123456789012345678901234567890123456789012


3962: c#######################################################################

3964:       implicit none

3966:       double precision fd, fc, fu

3968: c#######################################################################

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

3972:       return
3973:       end
3974:       double precision function  rexact(x,t)

3976:       implicit none

3978:       double precision x,t
3979:       double precision xot, head, tail, contact, ufan
3980:       double precision xpow, grat, urat
3981:       double precision uexact


3984:       logical debug

3986:       include 'tube.h'

3988:       debug = .false.


3991:       if (t .le. 0.0d+0) then
3992:         if (x .gt. 0.0d+0) then
3993:           rexact = r1
3994:         else
3995:           rexact = r4
3996:         endif
3997:       else

3999:        xot = x/t
4000:        head = -a4
4001:        tail = v3 - a3
4002:        contact = v2

4004:        if (xot .lt. head) then
4005:           rexact = r4
4006:        elseif (xot .gt. sspd) then
4007:           rexact = r1
4008:        elseif (xot .gt. contact) then
4009:           rexact = r2
4010:        elseif (xot .gt. tail) then
4011:           rexact = r3
4012:        else
4013:           ufan = uexact(x,t)
4014:           grat = (gamma - 1.0d+0) / 2.0d+0
4015:           xpow = 1.0d+0 / grat
4016:           urat = ufan / a4
4017:           rexact = r4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
4018:        endif

4020:       endif


4023:       if (debug) then
4024:         write(*,*)
4025:         write(*,*) 'rexact(',x,',',t,') = ',rexact
4026:         write(*,*)
4027:       endif

4029:       return
4030:       end
4031:       subroutine roestat(uhl, uhr, hl,hr,ht,uht)
4032: c23456789012345678901234567890123456789012345678901234567890123456789012
4033: c
4034: c          subroutine roestat
4035: c
4036: c  This subroutine computes the roe state at a face
4037: c
4038: c23456789012345678901234567890123456789012345678901234567890123456789012


4041: c#######################################################################

4043:       implicit none

4045:       include 'comd.h'

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

4049:       double precision ul, ur, shl, shr, xnum, xdenom
4050: 


4053: c#######################################################################

4055:       ul = uhl / hl
4056:       ur = uhr / hr

4058:       shl = sqrt(hl)
4059:       shr = sqrt(hr)

4061:       xnum = (shl * ul) + (shr * ur)
4062:       xdenom = shl + shr

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

4067:       return
4068:       end
4069:       subroutine rval2

4071:       implicit none

4073:       double precision prat, grat, xnum, xdenom


4076:       logical debug

4078:       include 'tube.h'

4080:       debug = .false.

4082:       prat = p2/p1
4083:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

4085:       xnum = 1.0d+0 + (grat * prat)
4086:       xdenom = grat + prat
4087: 
4088:       r2 = r1 * (xnum/xdenom)
4089: 


4092:       if (debug) then
4093:         write(*,*)
4094:         write(*,*) 'r1  = ',r1
4095:         write(*,*) 'r2  = ',r2
4096:       endif

4098:       return
4099:       end
4100:       subroutine  secondq
4101:      &      (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
4102:      &               rrg, rlg,rurg, rulg, erg, elg)
4103: c23456789012345678901234567890123456789012345678901234567890123456789012
4104: c
4105: c          subroutine secondq
4106: c
4107: c  this subroutine computes the second order (based on quick) left
4108: c  and right states for the godunov solver.
4109: c
4110: c23456789012345678901234567890123456789012345678901234567890123456789012


4113: c#######################################################################

4115:       implicit none

4117:       include 'comd.h'

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



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

4126:       double precision fluxlim2


4129: c#######################################################################

4131: c
4132: c  compute the velocities
4133: c
4134:       ull = rull/rll
4135:       ul  = rul /rl
4136:       ur  = rur /rr
4137:       urr = rurr/rrr

4139: c
4140: c  compute the left state first
4141: c
4142:       veld = 1.0d+0

4144:       rlg = fluxlim2(rll,rl,rr,rrr,veld)
4145:       ulg = fluxlim2(ull,ul,ur,urr,veld)
4146:       rulg = rlg * ulg
4147:       elg = fluxlim2(ell,el,er,err,veld)
4148: c
4149: c  now compute the right state
4150: c
4151:       veld = -1.0d+0

4153:       rrg = fluxlim2(rll,rl,rr,rrr,veld)
4154:       urg = fluxlim2(ull,ul,ur,urr,veld)
4155:       rurg = rrg * urg
4156:       erg = fluxlim2(ell,el,er,err,veld)



4160:       return
4161:       end
4162:       double precision function shockp(x)

4164:       implicit none

4166:       double precision x
4167:       double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
4168:       logical debug

4170:       include 'tube.h'

4172:       debug = .false.


4175:       if (debug) then
4176:          write(*,*)
4177:          write(*,*) 'gamma = ',gamma
4178:          write(*,*) 'a1 = ',a1
4179:          write(*,*) 'a4 = ',a4
4180:          write(*,*) 'p1 = ',p1
4181:          write(*,*) 'p2 = ',x
4182:          write(*,*)
4183:       endif

4185:       xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
4186:       xdenom = sqrt  (  2.0d+0 * gamma * ( (2.0d+0*gamma)
4187:      &            + (gamma + 1.0d+0) * ((x/p1) - 1) )  )
4188:       xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)

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


4193:       if (debug) then
4194:          write(*,*)
4195:          write(*,*) 'xnum = ',xnum
4196:          write(*,*) 'gamma = ',gamma
4197:          write(*,*) 'a1 = ',a1
4198:          write(*,*) 'a4 = ',a4
4199:          write(*,*) 'p1 = ',p1
4200:          write(*,*) 'xdenom = ',xdenom
4201:          write(*,*) 'xpow = ',xpow
4202:          write(*,*) 'shockp = ',shockp
4203:          write(*,*) 'p2 = ',x
4204:          write(*,*)
4205:       endif

4207:       return
4208:       end
4209:       double precision function  uexact(x,t)

4211:       implicit none

4213:       double precision x,t
4214:       double precision xot, head, tail


4217:       logical debug

4219:       include 'tube.h'

4221:       debug = .false.

4223:       if (debug) then
4224:         write(*,*)
4225:         write(*,*) 't = ',t
4226:         write(*,*) 'x = ',x
4227:         write(*,*) 'a4 = ',a4
4228:         write(*,*) 'v3 = ',v3
4229:         write(*,*) 'a3 = ',a3
4230:         write(*,*)
4231:       endif

4233:       if (t .le. 0.0d+0) then
4234:         uexact = 0.0d+0
4235:       else

4237:        xot = x/t
4238:        head = -a4
4239:        tail = v3 - a3

4241:        if (xot .lt. head) then
4242:           uexact = 0.0d+0
4243:        elseif (xot .gt. sspd) then
4244:           uexact = 0.0d+0
4245:        elseif (xot .gt. tail) then
4246:           uexact = v2
4247:        else
4248:           uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
4249:        endif

4251:       endif


4254:       if (debug) then
4255:         write(*,*)
4256: CVAM    write(*,*) 'x = ',x,' t = ',t
4257:         write(*,*) 'uexact = ',uexact
4258:         write(*,*)
4259:       endif

4261:       return
4262:       end
4263:       double precision function upwind(fw, fe, vp)
4264: c23456789012345678901234567890123456789012345678901234567890123456789012
4265: c
4266: c          function upwind
4267: c
4268: c  this function computes the upwind face value
4269: c
4270: c23456789012345678901234567890123456789012345678901234567890123456789012


4273: c#######################################################################

4275:       implicit none

4277:       double precision fw, fe, vp

4279: c#######################################################################

4281:       if (vp .gt. 0.0) then
4282:          upwind = vp * fw
4283:       else
4284:          upwind = vp * fe
4285:       endif

4287:       return
4288:       end
4289:       subroutine uval2

4291:       implicit none

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


4296:       logical debug

4298:       include 'tube.h'

4300:       debug = .false.

4302:       prat = p2/p1
4303:       grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
4304:       grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
4305:       arat = a1/gamma

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

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

4311:       if (debug) then
4312:         write(*,*)
4313:         write(*,*) 'v2  = ',v2
4314:       endif

4316:       return
4317:       end
4318:       subroutine val3

4320:       implicit none

4322:       double precision prat, rpow, epow, p3t


4325:       logical debug

4327:       include 'tube.h'

4329:       debug = .false.


4332:       p3 = p2

4334:       prat = p3/p4

4336:       rpow = 1.0d+0 / gamma

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

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

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

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

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

4348:       if (debug) then
4349:         write(*,*)
4350:         write(*,*) 'a3 = ',a3
4351:         write(*,*) 'r3 = ',r3
4352:         write(*,*) 'e3 = ',e3
4353:         write(*,*) 'p3 = ',p3
4354:         write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
4355:         write(*,*)
4356:       endif

4358:       return
4359:       end
4360:       subroutine wval

4362:       implicit none

4364:       double precision prat, grat, xnum


4367:       logical debug

4369:       include 'tube.h'

4371:       debug = .false.

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

4376:       xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
4377: 
4378:       sspd = a1 * sqrt(xnum)
4379: 


4382:       if (debug) then
4383:         write(*,*)
4384:         write(*,*) 'sspd  = ',sspd
4385:       endif

4387:       return
4388:       end