Actual source code: ex74f.F90

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

 37:       implicit none

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

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

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

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


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

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

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

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


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

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

108: !  Initialize problem parameters

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

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

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

134:       zero = 0.0
135:       ngraph = 10

137:       ndb = 0
138:       npc = 1

140:       damfac = 0.9d+0

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

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

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

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

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

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

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

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

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

171: !
172: !  load the command line options
173: !
174:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
175:      &                         '-dt',dt,flg,ierr)
176:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
177:      &                        '-mx',mx,flg,ierr)
178:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
179:      &                        '-nstep',nstep,flg,ierr)
180:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
181:      &                        '-debug',ndb,flg,ierr)
182:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
183:      &                        '-pcnew',npc,flg,ierr)
184:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
185:      &                        '-ihod',ihod,flg,ierr)
186:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
187:      &                        '-ientro',ientro,flg,ierr)
188:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
189:      &                         '-theta',theta,flg,ierr)
190:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
191:      &                        '-ngraph',ngraph,flg,ierr)
192:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
193:      &                         '-damfac',damfac,flg,ierr)
194:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
195:      &                        '-dampit',ndamp,flg,ierr)
196:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
197:      &                        '-wilson',nwilson,flg,ierr)
198:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
199:      &                        '-gorder',gorder,flg,ierr)
200:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
201:      &                        '-probnum',probnum,flg,ierr)
202:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
203:      &                         '-kappa0',kappa0,flg,ierr)
204:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
205:      &                         '-erg0',erg0,flg,ierr)
206:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
207:      &                         '-dtcon',ndtcon,flg,ierr)
208:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
209:      &                         '-tfinal',tfinal,flg,ierr)
210:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
211:      &                         '-tplot',tplot,flg,ierr)
212:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
213:      &                         '-dtgrow',dtgrow,flg,ierr)
214:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
215:      &                         '-tcscal',tcscal,flg,ierr)
216:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
217:      &                         '-hcscal',hcscal,flg,ierr)
218:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
219:      &                         '-dtmax',dtmax,flg,ierr)

221:       if (ndamp .eq. 1) then
222:          dampit = .true.
223:       endif

225:       if (nwilson .eq. 0) then
226:          wilson = .false.
227:       endif

229:       if (ndb .eq. 1) then
230:          debug = .true.
231:       endif

233:       if (npc .eq. 0) then
234:          pcnew = .false.
235:       endif

237:       if (ndtcon .eq. 0) then
238:          dtcon = .false.
239:       endif

241: !CVAM  if (dt .ge. dtmax .or. dt .le. dtmin) then
242: !CVAM     if (rank .eq. 0) write(6,*) 'DT is out of range'
243: !CVAM     SETERRA(1,0,' ')
244: !CVAM  endif

246:       N       = mx*neq

248:       ctx(5) = mx
249:       ctx(6) = N

251:       if (debug) then
252:         write(*,*) 'mx = ',mx
253:       endif



257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: !  Create nonlinear solver context
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

261:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: !  Create vector data structures; set function evaluation routine
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

269:       if (debug) then
270:         call MatGetSize(ctx(2),nx,ny,ierr)
271:         write(*,*) 'number of rows = ',nx,' number of col = ',ny
272:       endif
273: !
274: !  full size vectors
275: !
276:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,x,ierr)
277:       call VecSetFromOptions(x,ierr)
278:       call VecDuplicate(x,r,ierr)
279:       call VecDuplicate(x,ctx(4),ierr)
280: !
281: ! set grid
282: !
283:       dx = 2.0d+0/dfloat(mx)
284:       xl0 = -1.0d+0 -(0.5d+0 * dx)

286:       if (debug) then
287:         write(*,*) 'dx = ',dx
288:       endif
289: 

291: !  Set function evaluation routine and vector.  Whenever the nonlinear
292: !  solver needs to evaluate the nonlinear function, it will call this
293: !  routine.
294: !   - Note that the final routine argument is the user-defined
295: !     context that provides application-specific data for the
296: !     function evaluation routine.

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

300: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301: !  Customize nonlinear solver; set runtime options
302: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

306:       call SNESSetFromOptions(snes,ierr)
307: !
308: !  set the line search function to damp the newton update.
309: !
310:       if (dampit) then
311:         call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
312:       endif
313: !
314: ! get the linear solver info from the nonlinear solver
315: !

317:       call SNESGetKSP(snes,ksp,ierr)
318:       call KSPGetPC(ksp,pc,ierr)

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

323:       write(*,*)
324:       write(*,*)
325:       write(*,*) 'Linear solver'
326:       write(*,*)
327:       write(*,*) 'rtol = ',krtol
328:       write(*,*) 'atol = ',katol
329:       write(*,*) 'dtol = ',kdtol
330:       write(*,*) 'maxits = ',kmit
331:       write(*,*)
332:       write(*,*)
333:       write(*,*) 'Nonlinear solver'
334:       write(*,*)
335:       write(*,*) 'rtol = ',nrtol
336:       write(*,*) 'atol = ',natol
337:       write(*,*) 'stol = ',nstol
338:       write(*,*) 'maxits = ',nmit
339:       write(*,*) 'max func = ',nmf
340:       write(*,*)
341:       write(*,*)

343: !
344: !  Build shell based preconditioner if flag set
345: !
346:       if (pcnew) then
347:         call PCSetType(pc,PCSHELL,ierr)
348:         call PCShellSetContext(pc,ctx,ierr)
349:         call PCShellSetSetUpCtx(pc,PCRadSetUp,ierr)
350:         call PCShellSetApplyCtx(pc,PCRadApply,ctx,ierr)
351:       endif

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

355: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: !  Evaluate initial guess; then solve nonlinear system.
357: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358: !
359: !  initial counters
360: !
361:       time = 0.0d+0
362:       plotim = 0.0d+0
363:       totits = 0
364:       totlits = 0

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

371:       call FormInitialGuess(x,ierr)
372: !
373: !  open a window to plot results
374: !
375:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'density',0,0,300,300,view0,ierr)
376:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'velocity',320,0,300,300,view1,ierr)
377:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'total energy',640,0,300,300,view2,ierr)
378:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'temperature',0,360,300,300,view3,ierr)
379:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'pressure',320,360,300,300,view4,ierr)
380: !
381: !  graph initial conditions
382: !
383:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)
384: !
385: !  copy x into xold
386: !
387:       call VecCopy(x,ctx(4),ierr)
388:       call FormDt(snes,x,ctx,ierr)
389: !
390: !################################
391: !
392: !  TIME STEP LOOP BEGIN
393: !
394: !################################
395: !
396:       ndt = 0

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

400:         if (debug) then
401:           write(*,*)
402:           write(*,*) 'start of time loop'
403:           write(*,*)
404:           write(*,*) 'ndt = ',ndt
405:           write(*,*) 'nstep = ',nstep
406:           write(*,*) 'time = ',time
407:           write(*,*) 'tfinal = ',tfinal
408:           write(*,*)
409:         endif

411:         ndt = ndt + 1
412: !
413: !  increment time
414: !
415:         time = time + dt
416:         plotim = plotim + dt
417: !
418: !  call the nonlinear solver
419: !
420:         call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
421: !
422: !  get the number of linear iterations used by the nonlinear solver
423: !
424:         call SNESGetLinearSolveIterations(snes,lits,ierr)
425:         call SNESGetIterationNumber(snes,its,ierr)

427:         if (debug) then
428:            write(*,*) 'in radhyd ',ndt,'x'
429:            call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
430:         endif
431: !
432: !  increment the counters
433: !
434:         totits = totits + its
435:         totlits = totlits + lits
436: !
437: !  Compute new time step
438: !
439:           call FormDt(snes,x,ctx,ierr)
440: !
441: !  Draw contour plot of solution
442: !
443:         if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
444: 
445:            plotim = plotim - tplot


448:         if (rank .eq. 0) then
449:            write(6,100) totits,totlits,ndt,dt,time
450:         endif
451:   100   format('Newt = ',i7,' lin =',i7,' step =',i7,' dt = ',e8.3,' time = ',e10.4)
452: !
453: !  graph state conditions
454: !
455:           call FormGraph(x,view0,view1,view2,view3,view4,ierr)

457:         endif
458: !
459: ! copy x into xold
460: !
461:         call VecCopy(x,ctx(4),ierr)


464:         goto 10

466:       endif
467: !
468: !################################
469: !
470: !  TIME STEP LOOP END
471: !
472: !################################
473: !

475: !
476: !  graph final conditions
477: !
478:       call FormGraph(x,view0,view1,view2,view3,view4,ierr)


481:       write(*,*)
482:       write(*,*)
483:       write(*,*) 'total Newton iterations = ',totits
484:       write(*,*) 'total linear iterations = ',totlits
485:       write(*,*) 'Average Newton per time step = ', dble(totits)/dble(ndt)
486:       write(*,*) 'Average Krylov per Newton = ', dble(totlits)/dble(totits)
487:       write(*,*)
488:       write(*,*)

490: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
491: !  Free work space.  All PETSc objects should be destroyed when they
492: !  are no longer needed.
493: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


496:       call MatDestroy(ctx(2),ierr)
497:       call VecDestroy(x,ierr)
498:       call VecDestroy(ctx(4),ierr)
499:       call VecDestroy(r,ierr)
500:       call SNESDestroy(snes,ierr)
501:       call PetscViewerDestroy(view0,ierr)
502:       call PetscViewerDestroy(view1,ierr)
503:       call PetscViewerDestroy(view2,ierr)
504:       call PetscViewerDestroy(view3,ierr)
505:       call PetscViewerDestroy(view4,ierr)

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

509:       call PetscFinalize(ierr)

511:       close(87)

513:       stop
514:       end
515:       subroutine ApplicationDampit(x,deltx,w,ierr)
516: ! ---------------------------------------------------------------------
517: !
518: !  ApplicationDampit - Damps the newton update, called by
519: !  the higher level routine FormDampit().
520: !
521: !  Input Parameter:
522: !  x    - current iterate
523: !  deltx   - update
524: !
525: !  Output Parameters:
526: !  x    - new iterate
527: !  ierr - error code
528: !
529: !  Notes:
530: !  This routine only damps the density.  May want to add energy
531: !  in the future
532: !

534:       implicit none

536: !  Common blocks:
537: #include "ex74fcomd.h"

539: !  Input/output variables:
540:       PetscScalar   x(mx*neq), deltx(mx*neq), w(mx*neq)
541:       integer  ierr

543: !  Local variables:
544:       double precision facmin, fac, newx, xmin, se, dse
545:       double precision u,en,rn,run
546:       integer  i, jr, jru, je

548:       0

550:       if (debug) then
551:         write(*,*) 'begin damping'
552:         do i = 1,mx*neq
553:           write(*,*)'x(',i,') = ',x(i)
554:         enddo
555:         write(*,*)
556:         do i = 1,mx*neq
557:           write(*,*)'deltx(',i,') = ',deltx(i)
558:         enddo
559:       endif

561:       facmin = 1.0d+0
562: !
563: !  set the scale factor
564: !
565:       do i=1,mx
566: !
567: !  set pointers
568: !
569:         jr  = (neq*i) - 2
570:         jru = (neq*i) - 1
571:         je  = (neq*i)
572: !
573: !  make sure dencity stayes positive
574: !
575:         newx = x(jr) - deltx(jr)
576:         xmin = damfac * x(jr)

578:         if (newx .lt. xmin) then
579:           fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
580:           if (fac .lt. facmin) then
581:             if (debug) then
582:               write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
583:             endif
584:             facmin = fac
585:           endif
586:         endif
587: !
588: !  make sure Total energy stayes positive
589: !
590:         newx = x(je) - deltx(je)
591:         xmin = damfac * x(je)

593:         if (newx .lt. xmin) then
594:           fac = (1.0d+0 - damfac)*x(je)/deltx(je)
595:           if (fac .lt. facmin) then
596:             if (debug) then
597:               write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
598:             endif
599:             facmin = fac
600:           endif
601:         endif
602: !
603: !  make sure specific internal  energy stayes positive
604: !
605: 
606:         u = x(jru)/x(jr)
607:         se = (x(je)/x(jr)) - (0.5d+0 * u * u)

609:         en = x(je) - deltx(je)
610:         rn = x(jr) - deltx(jr)
611:         run = x(jru) - deltx(jru)

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


616:         newx = se - dse
617:         xmin = damfac * se

619:         if (newx .lt. xmin) then
620:           fac = (1.0d+0 - damfac) * se / dse
621:           if (fac .lt. facmin) then
622:             if (debug) then
623:               write(*,*) 'se',i, damfac,facmin,fac,se,dse
624:             endif
625:             facmin = fac
626:           endif
627:         endif

629:       enddo
630: !
631: ! write out warning
632: !
633:       if (facmin .lt. 1.0d+0) then
634:         write(*,*) 'facmin = ',facmin, damfac,time
635: !
636: !  scale the vector
637: !
638:         do i=1,neq*mx
639:           w(i) = x(i) - (facmin * deltx(i))
640:         enddo
641:       else
642:         do i=1,neq*mx
643:           w(i) = x(i) -  deltx(i)
644:         enddo
645:       endif

647:       if (debug) then
648:         write(*,*) 'end damping'
649:         do i = 1,mx*neq
650:            write(*,*) 'w(',i,') = ',w(i)
651:         enddo
652:       endif

654:       return
655:       end
656:       subroutine ApplicationDt(x,xold,ierr)
657: ! ---------------------------------------------------------------------
658: !
659: !  ApplicationDt - compute CFL numbers. Called by
660: !  the higher level routine FormDt().
661: !
662: !  Input Parameter:
663: !  x    - local vector data
664: !
665: !  Output Parameters:
666: !  ierr - error code
667: !
668: !  Notes:
669: !  This routine uses standard Fortran-style computations over a 2-dim array.
670: !

672:       implicit none

674: !  Common blocks:
675: #include "ex74fcomd.h"
676: #include "ex74ftube.h"

678: !  Input/output variables:
679:       PetscScalar   x(mx*neq), xold(mx*neq)
680:       integer  ierr

682: !  Local variables:
683:       integer  i, jr, jru, je
684: !
685: ! new
686: !
687:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw
688:       double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
689:       double precision tcflg, hcflg, dtt, dth
690:       double precision te, tp, tw
691:       double precision ue, up, uw
692:       double precision see, sep, sew
693: !
694: ! old
695: !
696:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
697:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
698:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
699:       double precision veloe,  velop,  velow
700:       double precision uop, seop, top
701:       double precision dtold, dttype
702: !
703: !  functions
704: !
705:       double precision eos

707:       dtold = dt

709:       0

711:       if (debug) then
712:         write(*,*) 'in start dt'
713:         do i = 1,mx*neq
714:           write(*,*)'x(',i,') = ',x(i)
715:         enddo
716:         write(*,*) 'tfinal = ',tfinal
717:         write(*,*) 'time = ',time
718:         write(*,*) 'dt = ',dt
719:         write(*,*) 'dtmax = ',dtmax
720:       endif

722:       sndp = -1.0d+20
723:       vradn = 0.0d+0
724:       vradd = 0.0d+0

726: !
727: !################################
728: !
729: !  loop over all cells begin
730: !
731: !################################
732: !
733:       do i=1,mx
734: !
735: !  set pointers
736: !
737:         jr  = (neq*i) - 2
738:         jru = (neq*i) - 1
739:         je  = (neq*i)
740: !
741: !
742: !  set scalars
743: !
744:         call Setpbc(i,x, rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww,vele,  velp,  velw)
745: !
746: ! compute temperatures
747: !
748:         uw = rhouw / rhow
749:         up = rhoup / rhop
750:         ue = rhoue / rhoe

752:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
753:         sep = (ergp/rhop) - (0.5d+0 * up * up)
754:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

756:         te  = see / csubv
757:         tp  = sep / csubv
758:         tw  = sew / csubv
759: !
760: ! compute old temperature
761: !

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

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

767:         uop = rhouop / rhoop

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

771:         top  = seop / csubv
772: !
773: !  compute thermal cfl
774: !
775:         vradn = vradn + (abs(tp - top)/dt)
776:         vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
777: !
778: !  compute hydro cfl
779: !

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

784:       enddo
785: !
786: !################################
787: !
788: !  loop over all cells end
789: !
790: !################################
791: !

793:       vrad = vradn / vradd

795:       tcfl = (vrad * dt) / dx
796:       hcfl = (sndp * dt) / dx

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

801:       dth = hcscal * dx / sndp

803:       if (.not. dtcon) then
804:         dt = min (dth,dtt,dt*dtgrow)
805:         if (dt .lt. dtmin) then
806:            dt = dtmin
807:         endif
808:         if (dt .gt. dtmax) then
809:            dt = dtmax
810:         endif
811:         if ( (time + dt) .gt. tfinal) then
812:            dt = tfinal - time
813:         endif

815:         if (dt .eq. dth) then
816:            dttype = 1.0d+0
817:         elseif (dt .eq. dtt) then
818:            dttype = 2.0d+0
819:         elseif (dt .eq. dtold*dtgrow) then
820:            dttype = 3.0d+0
821:         elseif (dt .eq. dtmax) then
822:            dttype = 4.0d+0
823:         elseif (dt .eq. dtmin) then
824:            dttype = 5.0d+0
825:         elseif (dt .eq. tfinal - time) then
826:            dttype = 6.0
827:         else
828:            dttype = -1.0d+0
829:         endif

831:       endif
832: 
833: 
834:       write(87,1000) time,dt,dth/hcscal,dtt/tcscal
835:       write(88,1000) time,dttype

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

839:       if (debug) then
840:         write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
841:         write(*,*) 'dtt = ',dtt,' dth = ',dth
842:         write(*,*) 'tfinal = ',tfinal
843:         write(*,*) 'time = ',time
844:         write(*,*) 'dt = ',dt
845:         write(*,*) 'dtmax = ',dtmax
846:         write(*,*)
847:       endif

849:       return
850:       end
851:       subroutine ApplicationExact(x,ierr)
852: ! ---------------------------------------------------------------------
853: !
854: !  ApplicationExact - Computes exact solution, called by
855: !  the higher level routine FormExact().
856: !
857: !  Input Parameter:
858: !  x - local vector data
859: !
860: !  Output Parameters:
861: !  x -    initial conditions
862: !  ierr - error code
863: !
864: !  Notes:
865: !  This routine uses standard Fortran-style computations over a 1-dim array.
866: !

868:       implicit none

870: !  Common blocks:

872: #include "ex74fcomd.h"

874: !  Input/output variables:
875:       PetscScalar  x(mx)
876:       integer ierr

878: !  Local variables:
879:       integer  i
880:       double precision xloc
881:       PetscScalar rexact


884: !  Set parameters

886:       0

888:       do i = 1,mx

890:         xloc = xl0 + (dble(i) * dx)
891:         x(i) = rexact(xloc,time)

893:       enddo

895:       return
896:       end
897:       subroutine ApplicationFunction(x,f,xold,ierr)
898: ! ---------------------------------------------------------------------
899: !
900: !  ApplicationFunction - Computes nonlinear function, called by
901: !  the higher level routine FormFunction().
902: !
903: !  Input Parameter:
904: !  x    - local vector data
905: !
906: !  Output Parameters:
907: !  f    - local vector data, f(x)
908: !  ierr - error code
909: !
910: !  Notes:
911: !  This routine uses standard Fortran-style computations over a 2-dim array.
912: !

914:       implicit none

916: !  Common blocks:
917: #include "ex74fcomd.h"

919: !  Input/output variables:
920:       PetscScalar   x(mx*neq), f(mx*neq), xold(mx*neq)
921:       integer  ierr

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

927:       double precision cont, energy, mom

929:       0

931:       if (debug) then
932:         write(*,*) 'in function'
933:         do i = 1,mx*neq
934:           write(*,*)'x(',i,') = ',x(i)
935:         enddo
936:       endif
937: !
938: !################################
939: !
940: !  loop over all cells begin
941: !
942: !################################
943: !
944:       do i=1,mx
945: !
946: !  set pointers
947: !
948:       jr  = (neq*i) - 2
949:       jru = (neq*i) - 1
950:       je  = (neq*i)
951: !
952: !
953: !  set scalars
954: !
955:       call Setpbc(i,x,rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, vele,  velp,  velw)
956: !
957: !  compute functions
958: !

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


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


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

968:        if (debug) then
969:          write(*,*)
970:          write(*,*) i,jr,jru,je,'res,r,ru,e'
971:          write(*,*) f(jr),f(jru),f(je)
972:          write(*,*)
973:        endif

975:       enddo
976: !
977: !################################
978: !
979: !  loop over all cells end
980: !
981: !################################
982: !

984:       if (debug) then
985:         write(*,*) 'in function'
986:         do i = 1,mx*neq
987:            write(*,*) 'f(',i,') = ',f(i)
988:         enddo
989:       endif

991:       return
992:       end
993:       subroutine ApplicationInitialGuess(x,ierr)
994: ! ---------------------------------------------------------------------
995: !
996: !  ApplicationInitialGuess - Computes initial approximation, called by
997: !  the higher level routine FormInitialGuess().
998: !
999: !  Input Parameter:
1000: !  x - local vector data
1001: !
1002: !  Output Parameters:
1003: !  x -    initial conditions
1004: !  ierr - error code
1005: !
1006: !  Notes:
1007: !  This routine uses standard Fortran-style computations over a 1-dim array.
1008: !

1010:       implicit none

1012: !  Common blocks:

1014: #include "ex74fcomd.h"
1015: #include "ex74ftube.h"

1017: !  Input/output variables:
1018:       PetscScalar  x(mx*neq)
1019:       integer ierr

1021: !  Local variables:
1022:       integer  i, j, jr, jru, je
1023:       double precision xloc, re, ee, ve
1024:       double precision wid, efloor
1025:       PetscScalar uexact, rexact, eexact


1028: !VAM  efloor = max(1.0d-1,1.0d-3 * erg0)
1029:       efloor = 1.0d-1
1030: !VAM  wid = max(1.0d-2,dx)
1031:       wid = 1.0d-2

1033: !  Set parameters

1035:       0

1037:       do i = 1,mx

1039:         jr  = (neq*i) - 2
1040:         jru = (neq*i) - 1
1041:         je  = (neq*i)

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

1045:         if (probnum .eq. 1) then
1046:            re = rexact(xloc,zero)
1047:            ve = uexact(xloc,zero)
1048:            ee = eexact(xloc,zero)
1049:         else
1050:            re = 1.0d+0
1051:            ve = 0.0d+0
1052:            ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1053:         endif

1055:         x(jr)  = re
1056:         x(jru) = re * ve
1057:         x(je)  = re * ( (0.5d+0 * ve * ve) + ee )

1059:         if (debug) then
1060:            write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1061:  100       format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1062:         endif

1064:       enddo

1066:       call exact0
1067:       call eval2
1068:       call rval2
1069:       call wval
1070:       call uval2
1071:       v3 = v2
1072:       call val3

1074:       a1 = sqrt(gamma*p1/r1)
1075:       a2 = sqrt(gamma*p2/r2)
1076:       a3 = sqrt(gamma*p3/r3)
1077:       a4 = sqrt(gamma*p4/r4)

1079:       write(*,1000) r1,r2,r3,r4
1080:       write(*,2000) p1,p2,p3,p4
1081:       write(*,3000) e1,e2,e3,e4
1082:       write(*,4000) a1,a2,a3,a4
1083:       write(*,*)

1085:  1000 format ('rhos      ',4(f12.6))
1086:  2000 format ('pressures ',4(f12.6))
1087:  3000 format ('energies  ',4(f12.6))
1088:  4000 format ('sound     ',4(f12.6))


1091:       return
1092:       end
1093:       subroutine ApplicationXmgr(x,ivar,ierr)
1094: ! ---------------------------------------------------------------------
1095: !
1096: !  ApplicationXmgr - Sets the Xmgr output called from
1097: !  the higher level routine FormXmgr().
1098: !
1099: !  Input Parameter:
1100: !  x - local vector data
1101: !
1102: !  Output Parameters:
1103: !  x -    initial conditions
1104: !  ierr - error code
1105: !
1106: !  Notes:
1107: !  This routine uses standard Fortran-style computations over a 1-dim array.
1108: !

1110:       implicit none

1112: !  Common blocks:

1114: #include "ex74fcomd.h"

1116: !  Input/output variables:
1117:       PetscScalar  x(mx)
1118:       integer ivar,ierr

1120: !  Local variables:
1121:       integer  i
1122:       double precision xloc, sum
1123:       PetscScalar rexact
1124:       integer iplotnum(5)
1125:       save iplotnum
1126:       character*8 grfile


1129:       data iplotnum / -1,-1,-1,-1,-1 /



1133: !  Set parameters

1135:       iplotnum(ivar) = iplotnum(ivar) + 1
1136:       0

1138:       if (ivar .eq. 1) then
1139:          write(grfile,4000) iplotnum(ivar)
1140:  4000    format('Xmgrr',i3.3)
1141:       elseif (ivar .eq. 2) then
1142:          write(grfile,5000) iplotnum(ivar)
1143:  5000    format('Xmgru',i3.3)
1144:       elseif (ivar .eq. 3) then
1145:          write(grfile,6000) iplotnum(ivar)
1146:  6000    format('Xmgre',i3.3)
1147:       elseif (ivar .eq. 4) then
1148:          write(grfile,7000) iplotnum(ivar)
1149:  7000    format('Xmgrt',i3.3)
1150:       else
1151:          write(grfile,8000) iplotnum(ivar)
1152:  8000    format('Xmgrp',i3.3)
1153:       endif

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

1157:       do i = 1,mx

1159:         xloc = xl0 + (dble(i) * dx)
1160:         if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1161:           write(44,1000) xloc, x(i), rexact(xloc,time)
1162:         else
1163:           write(44,1000) xloc, x(i)
1164:         endif

1166:       enddo

1168:  1000 format(3(e18.12,2x))
1169:       close(44)

1171:       if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1172:         sum = 0.0d+0
1173:         do i = 1,mx
1174:            xloc = xl0 + (dble(i) * dx)
1175:            sum = sum + (x(i) - rexact(xloc,time)) ** 2
1176:         enddo
1177:         sum = sqrt(sum)

1179:         write(*,*)
1180:         write(*,*)  'l2norm of the density error is',sum
1181:         write(*,*)
1182:       endif


1185:       return
1186:       end
1187:       subroutine FormDampit(snes,ctx,x,f,g,y,w, fnorm,ynorm,gnorm,flag,ierr)
1188: ! ---------------------------------------------------------------------
1189: !
1190: !  FormDampit - damps the Newton update
1191: !
1192: !  Input Parameters:
1193: !  snes  - the SNES context
1194: !  x     - current iterate
1195: !  f     - residual evaluated at x
1196: !  y     - search direction (containes new iterate on output)
1197: !  w     - work vector
1198: !  fnorm - 2-norm of f
1199: !
1200: !  In this example the application context is a Fortran integer array:
1201: !      ctx(1) = shell preconditioner pressure matrix contex
1202: !      ctx(2) = semi implicit pressure matrix
1203: !      ctx(4) = xold  - old time values need for time advancement
1204: !      ctx(5) = mx    - number of control volumes
1205: !      ctx(6) = N     - total number of unknowns
1206: !
1207: !  Output Parameter:
1208: !  g     - residual evaluated at new iterate y
1209: !  y     - new iterate (contains search direction on input
1210: !  gnorm - 2-norm of g
1211: !  ynorm - 2-norm of search length
1212: !  flag  - set to 0 if the line search succeeds; -1 on failure
1213: !
1214: !  Notes:
1215: !  This routine serves as a wrapper for the lower-level routine
1216: !  "ApplicationDampit", where the actual computations are
1217: !  done using the standard Fortran style of treating the local
1218: !  vector data as a multidimensional array over the local mesh.
1219: !  This routine merely accesses the local vector data via
1220: !  VecGetArray() and VecRestoreArray().
1221: !
1222:       implicit none

1224:  #include petsc/finclude/petsc.h

1226: !  Input/output variables:
1227:       SNES             snes
1228:       Vec              x, f, g, y, w
1229:       PetscFortranAddr ctx(*)
1230:       PetscScalar           fnorm, ynorm, gnorm
1231:       integer          ierr, flag

1233: !  Common blocks:

1235: #include "ex74fcomd.h"

1237: !  Local variables:

1239: !  Declarations for use with local arrays:
1240:       PetscScalar      lx_v(0:1), ly_v(0:1), lw_v(0:1)
1241:       PetscOffset lx_i, ly_i , lw_i

1243: !
1244: !  set ynorm
1245: !
1246:       call VecNorm(y,NORM_2,ynorm,ierr)
1247: !
1248: !  copy x to w
1249: !
1250:       call VecCopy(x,w,ierr)
1251: !
1252: ! get pointers to x, y, w
1253: !

1255:       call VecGetArray(x,lx_v,lx_i,ierr)
1256:       call VecGetArray(y,ly_v,ly_i,ierr)
1257:       call VecGetArray(w,lw_v,lw_i,ierr)
1258: !
1259: !  Compute Damping
1260: !
1261:       call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1262: !
1263: !  Restore vectors x, y, w
1264: !
1265:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1266:       call VecRestoreArray(y,ly_v,ly_i,ierr)
1267:       call VecRestoreArray(w,lw_v,lw_i,ierr)
1268: !
1269: !  copy w to y
1270: !
1271:       call VecCopy(w,y,ierr)
1272: !
1273: !  compute new residual
1274: !
1275:       call SNESComputeFunction(snes,y,g,ierr)
1276:       call VecNorm(g,NORM_2,gnorm,ierr)
1277:       flag = 0

1279:       if (debug) then
1280:          write(*,*) 'form damp ynorm = ',ynorm
1281:          write(*,*) 'gnorm = ',gnorm
1282:       endif

1284:       return
1285:       end
1286:       subroutine FormDt(snes,x,ctx,ierr)
1287: ! ---------------------------------------------------------------------
1288: !
1289: !  FormDt - Compute CFL numbers
1290: !
1291: !  Input Parameters:
1292: !  snes  - the SNES context
1293: !  x     - input vector
1294: !
1295: !  In this example the application context is a Fortran integer array:
1296: !      ctx(1) = shell preconditioner pressure matrix contex
1297: !      ctx(2) = semi implicit pressure matrix
1298: !      ctx(4) = xold  - old time values need for time advancement
1299: !      ctx(5) = mx    - number of control volumes
1300: !      ctx(6) = N     - total number of unknowns
1301: !
1302: !
1303: !  Notes:
1304: !  This routine serves as a wrapper for the lower-level routine
1305: !  "ApplicationDt", where the actual computations are
1306: !  done using the standard Fortran style of treating the local
1307: !  vector data as a multidimensional array over the local mesh.
1308: !  This routine merely accesses the local vector data via
1309: !  VecGetArray() and VecRestoreArray().
1310: !
1311:       implicit none

1313:  #include petsc/finclude/petsc.h

1315: !  Input/output variables:
1316:       SNES             snes
1317:       Vec              x
1318:       PetscFortranAddr ctx(*)
1319:       integer          ierr

1321: !  Common blocks:

1323: #include "ex74fcomd.h"

1325: !  Local variables:

1327: !  Declarations for use with local arrays:
1328:       PetscScalar      lx_v(0:1)
1329:       PetscOffset lx_i
1330:       PetscScalar      lxold_v(0:1)
1331:       PetscOffset lxold_i

1333: !
1334: ! get pointers to x, xold
1335: !

1337:       call VecGetArray(x,lx_v,lx_i,ierr)
1338:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1339: !
1340: !  Compute function
1341: !
1342:       call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1343: !
1344: !  Restore vectors x, xold
1345: !
1346:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1347:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)

1349:       return
1350:       end
1351:       subroutine FormExact(x,ierr)
1352: ! ---------------------------------------------------------------------
1353: !
1354: !  FormExact - Forms exact solution
1355: !
1356: !  Input Parameter:
1357: !  x - vector
1358: !
1359: !  Output Parameters:
1360: !  x - vector
1361: !  ierr - error code
1362: !
1363: !  Notes:
1364: !  This routine serves as a wrapper for the lower-level routine
1365: !  "ApplicationExact", where the actual computations are
1366: !  done using the standard Fortran style of treating the local
1367: !  vector data as a multidimensional array over the local mesh.
1368: !  This routine merely accesses the local vector data via
1369: !  VecGetArray() and VecRestoreArray().
1370: !
1371:       implicit none

1373:  #include petsc/finclude/petsc.h

1375: !  Input/output variables:
1376:       Vec      x
1377:       integer  ierr

1379: !  Declarations for use with local arrays:
1380:       PetscScalar      lx_v(0:1)
1381:       PetscOffset lx_i

1383:       0

1385: !
1386: !  get a pointer to x
1387: !
1388:       call VecGetArray(x,lx_v,lx_i,ierr)
1389: !
1390: !  Compute initial guess
1391: !
1392:       call ApplicationExact(lx_v(lx_i),ierr)
1393: !
1394: !  Restore vector x
1395: !
1396:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1398:       return
1399:       end
1400:       subroutine FormFunction(snes,x,f,ctx,ierr)
1401: ! ---------------------------------------------------------------------
1402: !
1403: !  FormFunction - Evaluates nonlinear function, f(x).
1404: !
1405: !  Input Parameters:
1406: !  snes  - the SNES context
1407: !  x     - input vector
1408: !
1409: !  In this example the application context is a Fortran integer array:
1410: !      ctx(1) = shell preconditioner pressure matrix contex
1411: !      ctx(2) = semi implicit pressure matrix
1412: !      ctx(4) = xold  - old time values need for time advancement
1413: !      ctx(5) = mx    - number of control volumes
1414: !      ctx(6) = N     - total number of unknowns
1415: !
1416: !  Output Parameter:
1417: !  f     - vector with newly computed function
1418: !
1419: !  Notes:
1420: !  This routine serves as a wrapper for the lower-level routine
1421: !  "ApplicationFunction", where the actual computations are
1422: !  done using the standard Fortran style of treating the local
1423: !  vector data as a multidimensional array over the local mesh.
1424: !  This routine merely accesses the local vector data via
1425: !  VecGetArray() and VecRestoreArray().
1426: !
1427:       implicit none

1429:  #include petsc/finclude/petsc.h

1431: !  Input/output variables:
1432:       SNES             snes
1433:       Vec              x, f
1434:       PetscFortranAddr ctx(*)
1435:       integer          ierr

1437: !  Common blocks:

1439: #include "ex74fcomd.h"

1441: !  Local variables:

1443: !  Declarations for use with local arrays:
1444:       PetscScalar      lx_v(0:1), lf_v(0:1)
1445:       PetscOffset lx_i, lf_i
1446:       PetscScalar      lxold_v(0:1)
1447:       PetscOffset lxold_i

1449: !
1450: ! get pointers to x, f, xold
1451: !

1453:       call VecGetArray(x,lx_v,lx_i,ierr)
1454:       call VecGetArray(f,lf_v,lf_i,ierr)
1455:       call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1456: !
1457: !  Compute function
1458: !
1459:       call ApplicationFunction(lx_v(lx_i),lf_v(lf_i), lxold_v(lxold_i),ierr)
1460: !
1461: !  Restore vectors x, f, xold
1462: !
1463:       call VecRestoreArray(x,lx_v,lx_i,ierr)
1464:       call VecRestoreArray(f,lf_v,lf_i,ierr)
1465:       call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1466: !
1467: ! something to do with profiling
1468: !
1469:       call PetscLogFlops(110.d0*mx,ierr)

1471:       return
1472:       end
1473:       subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1474: ! ---------------------------------------------------------------------
1475: !
1476: !  FormGraph - Forms Graphics output
1477: !
1478: !  Input Parameter:
1479: !  x - vector
1480: !  time - scalar
1481: !
1482: !  Output Parameters:
1483: !  ierr - error code
1484: !
1485: !  Notes:
1486: !  This routine serves as a wrapper for the lower-level routine
1487: !  "ApplicationXmgr", where the actual computations are
1488: !  done using the standard Fortran style of treating the local
1489: !  vector data as a multidimensional array over the local mesh.
1490: !  This routine merely accesses the local vector data via
1491: !  VecGetArray() and VecRestoreArray().
1492: !
1493:       implicit none

1495:  #include petsc/finclude/petsc.h

1497: #include "ex74fcomd.h"
1498: #include "ex74ftube.h"

1500: !  Input/output variables:
1501:       Vec      x
1502:       integer  ierr

1504: !  Declarations for use with local arrays:
1505:       IS                 rfrom, rto, rufrom, ruto, efrom, eto
1506:       Vec                rval
1507:       Vec                uval
1508:       Vec                ruval
1509:       Vec                eval
1510:       Vec                seval
1511:       Vec                pval
1512:       Vec                kval
1513:       Vec                tval
1514:       Vec                steval
1515:       VecScatter         scatter
1516:       PetscViewer             view0, view1, view2, view3, view4
1517:       double precision minus1, l2err, gm1, csubvi


1520:       csubvi = 1.0d+0 / csubv
1521:       gm1 = gamma - 1.0d+0
1522:       0
1523:       minus1 = -1.0d+0
1524: !
1525: !  graphics vectors
1526: !
1527:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,rval,ierr)
1528:       call VecSetFromOptions(rval,ierr)
1529:       call VecDuplicate(rval,uval,ierr)
1530:       call VecDuplicate(rval,ruval,ierr)
1531:       call VecDuplicate(rval,eval,ierr)
1532:       call VecDuplicate(rval,seval,ierr)
1533:       call VecDuplicate(rval,pval,ierr)
1534:       call VecDuplicate(rval,kval,ierr)
1535:       call VecDuplicate(rval,tval,ierr)
1536:       call VecDuplicate(rval,steval,ierr)
1537: !
1538: ! create index sets for vector scatters
1539: !
1540:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1541:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  rto,   ierr)
1542:       call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1543:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  ruto,  ierr)
1544:       call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1545:       call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,  eto,   ierr)

1547: !
1548: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1549: !
1550: !
1551: !  load rval from x
1552: !
1553:       call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1554:       call VecScatterBegin(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1555:       call VecScatterEnd(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1556:       call VecScatterDestroy(scatter,ierr)
1557: !
1558: !  plot rval vector
1559: !
1560:       call VecView(rval,view0,ierr)
1561: !
1562: !  make xmgr plot of rval
1563: !
1564:       call FormXmgr(rval,1,ierr)
1565: !
1566: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1567: !
1568: !
1569: !  load eval from x
1570: !
1571:       call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1572:       call VecScatterBegin(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1573:       call VecScatterEnd(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1574:       call VecScatterDestroy(scatter,ierr)
1575: !
1576: !  plot eval vector
1577: !
1578:       call VecView(eval,view2,ierr)
1579: !
1580: !  make xmgr plot of eval
1581: !
1582:       call FormXmgr(eval,3,ierr)
1583: !
1584: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1585: !

1587: !
1588: !  load ruval from x
1589: !
1590:       call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1591:       call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1592:       call VecScatterEnd(scatter,x,ruval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1593:       call VecScatterDestroy(scatter,ierr)
1594: !
1595: !  create u = ru / r
1596: !
1597:       call VecPointwiseDivide(uval,ruval,rval,ierr)
1598: !
1599: !  plot uval vector
1600: !
1601:       call VecView(uval,view1,ierr)
1602: !
1603: !  make xmgr plot of uval
1604: !
1605:       call FormXmgr(uval,2,ierr)

1607: !
1608: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1609: !

1611:       call VecPointwiseMult(kval,uval,uval,ierr)
1612:       call VecScale(kval,0.5d+0,ierr)

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

1617:       call VecCopy(seval,tval,ierr)
1618:       call VecScale(tval,csubvi,ierr)

1620: !
1621: !  plot tval vector
1622: !
1623:       call VecView(tval,view3,ierr)
1624: !
1625: !  make xmgr plot of tval
1626: !
1627:       call FormXmgr(tval,4,ierr)

1629: !
1630: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1631: !

1633:       call VecPointwiseMult(pval,rval,seval,ierr)
1634:       call VecScale(pval,gm1,ierr)
1635: !
1636: !  plot pval vector
1637: !
1638:       call VecView(pval,view4,ierr)
1639: !
1640: !  make xmgr plot of pval
1641: !
1642:       call FormXmgr(pval,5,ierr)
1643: !
1644: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1645: !





1651: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1652: !  Free work space.  All PETSc objects should be destroyed when they
1653: !  are no longer needed.
1654: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1656:       call VecDestroy(rval, ierr)
1657:       call VecDestroy(uval, ierr)
1658:       call VecDestroy(ruval,ierr)
1659:       call VecDestroy(eval, ierr)
1660:       call VecDestroy(seval, ierr)
1661:       call VecDestroy(pval, ierr)
1662:       call VecDestroy(kval, ierr)
1663:       call VecDestroy(tval, ierr)
1664:       call VecDestroy(steval, ierr)

1666:       call ISDestroy(rfrom, ierr)
1667:       call ISDestroy(rto,   ierr)

1669:       call ISDestroy(rufrom,ierr)
1670:       call ISDestroy(ruto,  ierr)

1672:       call ISDestroy(efrom, ierr)
1673:       call ISDestroy(eto,   ierr)


1676:       return
1677:       end
1678:       subroutine FormInitialGuess(x,ierr)
1679: ! ---------------------------------------------------------------------
1680: !
1681: !  FormInitialGuess - Forms initial approximation.
1682: !
1683: !  Input Parameter:
1684: !  x - vector
1685: !
1686: !  Output Parameters:
1687: !  x - vector
1688: !  ierr - error code
1689: !
1690: !  Notes:
1691: !  This routine serves as a wrapper for the lower-level routine
1692: !  "ApplicationInitialGuess", where the actual computations are
1693: !  done using the standard Fortran style of treating the local
1694: !  vector data as a multidimensional array over the local mesh.
1695: !  This routine merely accesses the local vector data via
1696: !  VecGetArray() and VecRestoreArray().
1697: !
1698:       implicit none

1700:  #include petsc/finclude/petsc.h

1702: !  Input/output variables:
1703:       Vec      x
1704:       integer  ierr

1706: !  Declarations for use with local arrays:
1707:       PetscScalar      lx_v(0:1)
1708:       PetscOffset lx_i

1710:       0

1712: !
1713: !  get a pointer to x
1714: !
1715:       call VecGetArray(x,lx_v,lx_i,ierr)
1716: !
1717: !  Compute initial guess
1718: !
1719:       call ApplicationInitialGuess(lx_v(lx_i),ierr)
1720: !
1721: !  Restore vector x
1722: !
1723:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1725:       return
1726:       end
1727:       subroutine FormXmgr(x,ivar,ierr)
1728: ! ---------------------------------------------------------------------
1729: !
1730: !  FormXmgr - Forms Xmgr output
1731: !
1732: !  Input Parameter:
1733: !  x - vector
1734: !
1735: !  Output Parameters:
1736: !  x - vector
1737: !  ierr - error code
1738: !
1739: !  Notes:
1740: !  This routine serves as a wrapper for the lower-level routine
1741: !  "ApplicationXmgr", where the actual computations are
1742: !  done using the standard Fortran style of treating the local
1743: !  vector data as a multidimensional array over the local mesh.
1744: !  This routine merely accesses the local vector data via
1745: !  VecGetArray() and VecRestoreArray().
1746: !
1747:       implicit none

1749:  #include petsc/finclude/petsc.h

1751: !  Input/output variables:
1752:       Vec      x
1753:       integer  ivar,ierr

1755: !  Declarations for use with local arrays:
1756:       PetscScalar      lx_v(0:1)
1757:       PetscOffset lx_i

1759:       0

1761: !
1762: !  get a pointer to x
1763: !
1764:       call VecGetArray(x,lx_v,lx_i,ierr)
1765: !
1766: !  make the graph
1767: !
1768:       call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1769: !
1770: !  Restore vector x
1771: !
1772:       call VecRestoreArray(x,lx_v,lx_i,ierr)

1774:       return
1775:       end
1776:       subroutine PCRadApply(pc,ctx,x,y,ierr)
1777: ! -------------------------------------------------------------------
1778: !
1779: !   PCRadApply - This routine demonstrates the use of a
1780: !   user-provided preconditioner.
1781: !
1782: !   Input Parameters:
1783: !   dummy - optional user-defined context, not used here
1784: !   x - input vector
1785: !  In this example the shell preconditioner application context
1786: !  is a Fortran integer array:
1787: !      ctx(1) = shell preconditioner pressure matrix contex
1788: !      ctx(2) = semi implicit pressure matrix
1789: !      ctx(4) = xold  - old time values need for time advancement
1790: !      ctx(5) = mx    - number of control volumes
1791: !      ctx(6) = N     - total number of unknowns
1792: !
1793: !   Output Parameters:
1794: !   y - preconditioned vector
1795: !   ierr  - error code (nonzero if error has been detected)
1796: !
1797: !   Notes:
1798: !   This code implements the Jacobi preconditioner plus the
1799: !   SOR preconditioner
1800: !

1802:       implicit none

1804:  #include petsc/finclude/petsc.h

1806: #include "ex74fcomd.h"
1807: !
1808: !  Input
1809: !
1810:       PC               pc
1811:       PetscFortranAddr ctx(*)
1812:       Vec              x, y
1813:       integer          ierr
1814: !
1815: !  Local
1816: !
1817:       IS               defrom, deto
1818:       Vec              de, rese
1819:       VecScatter       scatter
1820:       PetscScalar           lde_v(0:1),lrese_v(0:1)
1821:       PetscOffset      lde_i,     lrese_i
1822: !
1823: !  Identity preconditioner
1824: !
1825:       call VecCopy(x,y,ierr)
1826: !
1827: !  if kappa0 not equal to zero then precondition the radiation diffusion
1828: !
1829:       if (kappa0 .ne. 0.0d+0) then
1830: 

1832: !
1833: !  Create needed vectors
1834: !
1835:          call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,de,ierr)
1836:          call VecSetFromOptions(de,ierr)
1837:          call VecDuplicate(de,rese,ierr)
1838: !
1839: !  create index sets for scatters
1840: !
1841:          call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1842:          call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1843: !
1844: !  load rese from x
1845: !
1846:          call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1847:          call VecScatterBegin(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1848:          call VecScatterEnd(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1849:          call VecScatterDestroy(scatter,ierr)
1850: !
1851: !  apply preconditioner
1852: !
1853:       call PCApply(ctx(1),rese,de,ierr)

1855:       if (debug) then
1856:         write(*,*) 'PCRadApply dh is'
1857:         call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1858:       endif
1859: !
1860: ! load de into y
1861: !
1862:       call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1863:       call VecScatterBegin(scatter,de,y,INSERT_VALUES, SCATTER_FORWARD,ierr)
1864:       call VecScatterEnd(scatter,de,y,INSERT_VALUES,SCATTER_FORWARD,ierr)
1865:       call VecScatterDestroy(scatter,ierr)

1867:       if (debug) then
1868:         write(*,*) 'PCRadApply y is'
1869:         call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1870:       endif

1872:       call VecDestroy(de,ierr)
1873:       call VecDestroy(rese,ierr)

1875:       call ISDestroy(defrom,ierr)
1876:       call ISDestroy(deto,ierr)

1878:       endif


1881:       return
1882:       end
1883:       subroutine PCRadSetUp(pc,ctx,ierr)
1884: !
1885: !   PCRadSetUp - This routine sets up a user-defined
1886: !   preconditioner context.
1887: !
1888: !   Input Parameters:
1889: !  In this example the shell preconditioner application context
1890: !  is a Fortran integer array:
1891: !      ctx(1) = shell preconditioner pressure matrix contex
1892: !      ctx(2) = semi implicit pressure matrix
1893: !      ctx(4) = xold  - old time values need for time advancement
1894: !      ctx(5) = mx    - number of control volumes
1895: !      ctx(6) = N     - total number of unknowns
1896: !
1897: !   Output Parameter:
1898: !   ierr  - error code (nonzero if error has been detected)
1899: !
1900: !   Notes:
1901: !   In this example, we define the shell preconditioner to be Jacobi
1902: !   method.  Thus, here we create a work vector for storing the reciprocal
1903: !   of the diagonal of the preconditioner matrix; this vector is then
1904: !   used within the routine PCRadApply().
1905: !

1907:       implicit none

1909:  #include petsc/finclude/petsc.h

1911: #include "ex74fcomd.h"
1912: !
1913: !  Input
1914: !
1915:       PC               pc
1916:       PetscFortranAddr ctx(*)
1917:       integer          ierr
1918: !
1919: !  Local
1920: !
1921:       Vec              eold
1922: 
1923:       PetscScalar      le_v(0:1)
1924:       PetscOffset le_i
1925: 
1926: !
1927: !  create vector
1928: !
1929:       call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,eold,ierr)
1930:       call VecSetFromOptions(eold,ierr)
1931: !
1932: ! set up the matrix based on xold
1933: !
1934:       call Setmat(ctx,ierr)
1935: !
1936: !  set up the preconditioner
1937: !
1938:       call PCDestroy(ctx(1),ierr)
1939:       call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
1940: !VAM  call PCSetType(ctx(1),PCJACOBI,ierr)
1941:       call PCSetType(ctx(1),PCLU,ierr)
1942: !      call PCSetVector(ctx(1),eold,ierr)
1943:       call PCSetOperators(ctx(1),ctx(2),ctx(2), ierr)
1944:       call PCSetUp(ctx(1),ierr)

1946:       call VecDestroy(eold,ierr)


1949:       return
1950:       end
1951:       subroutine Setmat(ctx,ierr)

1953:       implicit none

1955:  #include petsc/finclude/petsc.h

1957: !  Common blocks:
1958: #include "ex74fcomd.h"
1959: #include "ex74ftube.h"

1961: !  Input/output variables:
1962:       PetscFortranAddr ctx(*)
1963:       integer  ierr

1965: !  Local variables:
1966:       PetscScalar      lx_v(0:1)
1967:       PetscOffset lx_i

1969:       double precision xmult, himh, hiph, diag, upper, lower
1970:       double precision hi, hip1, him1
1971:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, ue,    up,    uw
1972:       double precision see, sep, sew, seef, sewf, tef, twf, ref, rwf, kef, kwf, xmulte, xmultw
1973: !
1974:       integer  im, nx, ny
1975: !
1976: !     get pointers to xold
1977: !
1978:       call VecGetArray(ctx(4),lx_v,lx_i,ierr)
1979: 

1981: !
1982: !############################
1983: !
1984: ! loop over all cells begin
1985: !
1986: !############################
1987: !
1988:       do im = 1,mx
1989: !
1990: !  set scalars
1991: !
1992:          call Setpbc(im,lx_v(lx_i), rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww, ue,    up,    uw)
1993: !
1994: !  set diffusion coefficients
1995: !
1996:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
1997:         sep = (ergp/rhop) - (0.5d+0 * up * up)
1998:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

2000:         seef = 0.5d+0 * (see + sep)
2001:         sewf = 0.5d+0 * (sew + sep)

2003:         tef = seef / csubv
2004:         twf = sewf / csubv

2006:         ref = 0.5d+0 * (rhoe + rhop)
2007:         rwf = 0.5d+0 * (rhow + rhop)

2009:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2010:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

2012:         if (wilson) then
2013:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2014:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2015:         endif
2016: !
2017: !  set coefficients
2018: !
2019:          xmult = dt / (dx * dx * csubv)

2021:          xmulte = xmult * kef
2022:          xmultw = xmult * kwf

2024:          upper = -(xmulte / rhoe)
2025:          lower = -(xmultw / rhow)

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

2029: !
2030: !  load coefficients into the matrix
2031: !
2032:          call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)

2034:          if (im .eq. 1) then
2035:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper, INSERT_VALUES,ierr)
2036:          elseif (im .eq. mx) then
2037:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower,INSERT_VALUES,ierr)
2038:          else
2039:            call MatSetValues(ctx(2),1,im-1,1,im  ,upper,INSERT_VALUES,ierr)
2040:            call MatSetValues(ctx(2),1,im-1,1,im-2,lower, INSERT_VALUES,ierr)
2041:          endif


2044:       enddo
2045: !
2046: !############################
2047: !
2048: ! loop over all cells end
2049: !
2050: !############################
2051: !
2052: 
2053: !
2054: !  final load of matrix
2055: !
2056:       call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2057:       call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)

2059:       if (debug) then
2060:         call MatGetSize(ctx(2),nx,ny,ierr)
2061:         write(*,*) 'in setup nx = ',nx,' ny = ',ny
2062:         call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2063:       endif

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



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

2073:       implicit none

2075: !  Common blocks:
2076: #include "ex74fcomd.h"

2078: !  Input/output variables:
2079:       PetscScalar   x(mx*neq)
2080:       integer  i
2081:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2082:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2083:       double precision ergee,  erge,  ergp,  ergw,  ergww
2084:       double precision         vele,  velp,  velw

2086: !  Local variables:
2087:       integer  jr, jru, je

2089: !
2090: !  set pointers
2091: !
2092:       jr  = (neq*i) - 2
2093:       jru = (neq*i) - 1
2094:       je  = (neq*i)

2096:       if (debug) then
2097:         write(*,*)
2098:         write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2099:         write(*,*)
2100:       endif
2101: 
2102:       if (i .eq. 1) then

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

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

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

2122:         velw = 0.0d+0
2123:         velp = rhoup/rhop
2124:         vele = rhoue/rhoe

2126:       elseif (i .eq. 2) then

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

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

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

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

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

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

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

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

2170:         velw = rhouw/rhow
2171:         velp = rhoup/rhop
2172:         vele = rhoue/rhoe

2174:       elseif (i .eq. mx) then

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

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

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

2194:         velw = rhouw/rhow
2195:         velp = rhoup/rhop
2196:         vele = 0.0d+0

2198:       else

2200:         rhoee = x(jr+(2*neq))
2201:         rhoe  = x(jr+neq)
2202:         rhop  = x(jr)
2203:         rhow  = x(jr-neq)
2204:         rhoww = x(jr-(2*neq))

2206:         rhouee = x(jru+(2*neq))
2207:         rhoue  = x(jru+neq)
2208:         rhoup  = x(jru)
2209:         rhouw  = x(jru-neq)
2210:         rhouww = x(jru-(2*neq))

2212:         ergee = x(je+(2*neq))
2213:         erge  = x(je+neq)
2214:         ergp  = x(je)
2215:         ergw  = x(je-neq)
2216:         ergww = x(je-(2*neq))

2218:         velw = rhouw/rhow
2219:         velp = rhoup/rhop
2220:         vele = rhoue/rhoe

2222:       endif

2224:       if (debug) then
2225:          write(*,*)
2226:          write(*,*) 'in Setpbc ',i,jr,jru,je
2227:          write(*,*) 'mx = ',mx
2228:          write(*,*)
2229:       endif


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

2236:       implicit none

2238: !  Common blocks:
2239: #include "ex74fcomd.h"

2241: !  Input/output variables:
2242:       integer  jbc
2243:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2244:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2245:       double precision ergee,  erge,  ergp,  ergw,  ergww
2246:       double precision         ue,    up,    uw

2248: !  Local variables:

2250:       if (jbc .eq. 1) then
2251:          rhoww  = rhop
2252:          rhow   = rhop
2253:          rhouww = rhoup
2254:          rhouw  = rhoup
2255:          ergww  = ergp
2256:          ergw   = ergp
2257:          uw     = 0.0d+0
2258:       elseif  (jbc .eq. 2) then
2259:          rhoww  = rhow
2260:          rhouww = rhouw
2261:          ergww  = ergw
2262:          uw     = rhouw / rhow
2263:       else
2264:          uw = rhouw / rhow
2265:       endif

2267:       if (jbc .eq. mx) then
2268:          rhoee  = rhop
2269:          rhoe   = rhop
2270:          rhouee = rhoup
2271:          rhoue  = rhoup
2272:          ergee  = ergp
2273:          erge   = ergp
2274:          ue     = 0.0d+0
2275:       elseif (jbc .eq. mx-1) then
2276:          rhoee  = rhoe
2277:          rhouee = rhoue
2278:          ergee  = erge
2279:          ue     = rhoue / rhoe
2280:       else
2281:          ue     = rhoue / rhoe
2282:       endif

2284:       up = rhoup / rhop

2286:       if (debug) then
2287:          write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2288:       endif


2291:       return
2292:       end
2293:       double precision function cont(rhoee,  rhoe,  rhop,  rhow,  rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee,  erge,  ergp,  ergw,  ergww,jcont,xold)
2294: !
2295: !  This function computes the residual
2296: !  for the 1-D continuity equation
2297: !
2298: !
2299:       implicit none

2301:       include 'ex74fcomd.h'
2302:       include 'ex74ftube.h'
2303: !
2304: !     input variables
2305: !
2306:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2307:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2308:       double precision ergee,  erge,  ergp,  ergw,  ergww
2309:       double precision xold(mx*neq)
2310: !
2311:       integer jcont
2312: !
2313: !     local variables
2314: !
2315:       double precision theta1
2316:       integer jr
2317: !
2318: !  new
2319: !
2320:       double precision velfw, velfe
2321:       double precision vele,velp,velw
2322:       double precision fluxe, fluxw
2323:       double precision urhoe, urhow
2324:       double precision source
2325: !
2326: ! old
2327: !
2328:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2329:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2330:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2331:       double precision teoee, teoe, teop, teow, teoww, uoe, uoee, uop, uow, uoww
2332:       double precision velfow, velfoe
2333:       double precision veloe,velop,velow
2334:       double precision fluxoe, fluxow
2335:       double precision urhooe, urhoow
2336:       double precision sourceo
2337: !
2338: ! functions
2339: !
2340:       double precision godunov2
2341:       double precision upwind, fluxlim
2342: !
2343: !
2344: ! ******************************************************************
2345: !
2346: !
2347: !
2348:       if (debug) then
2349:         write(*,*)
2350:         write(*,*) 'in cont',jcont,' ihod = ',ihod
2351:         write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2352:         write(*,*) 'rhop = ',rhop
2353:         write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2354:         write(*,*)
2355:       endif

2357:       jr = (neq*jcont) - 2

2359: !########################
2360: !
2361: !      NEW
2362: !
2363: !########################

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

2367:       velfe = 0.5d+0 * (vele + velp)
2368:       velfw = 0.5d+0 * (velw + velp)

2370:       if (ihod .eq. 1) then

2372:         urhoe = upwind(rhop,rhoe,velfe)
2373:         urhow = upwind(rhow,rhop,velfw)

2375:       elseif (ihod .eq. 2) then

2377:         urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2378:         urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)

2380:       endif

2382:       if (ihod .eq. 3) then
2383:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,1)
2384:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,1)
2385:       else
2386:         fluxe = (dt/dx) * urhoe
2387:         fluxw = (dt/dx) * urhow
2388:       endif


2391:       source = 0.0d+0

2393: !########################
2394: !
2395: !      OLD
2396: !
2397: !########################

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

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

2403:       velfoe = 0.5d+0 * (veloe + velop)
2404:       velfow = 0.5d+0 * (velow + velop)


2407:       if (ihod .eq. 1) then

2409:         urhooe = upwind(rhoop,rhooe,velfoe)
2410:         urhoow = upwind(rhoow,rhoop,velfow)

2412:       elseif (ihod .eq. 2) then

2414:         urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2415:         urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)

2417:       endif

2419:       if (ihod .eq. 3) then
2420:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,1)
2421:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,1)
2422:       else
2423:         fluxoe = (dt/dx) * urhooe
2424:         fluxow = (dt/dx) * urhoow
2425:       endif

2427:       sourceo = 0.0d+0


2430: !########################
2431: !
2432: !      FUNCTION
2433: !
2434: !########################

2436:       theta1 = 1.0d+0 - theta
2437:       cont =  (rhop - xold(jr))  + (  theta  * ( (fluxe  - fluxw ) - source  )  )   + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
2438: !VAM
2439:       if (probnum .eq. 3) then
2440:         cont = 0.0d+0
2441:       endif
2442: !VAM


2445:       if (debug) then
2446:        write(*,*)
2447:        write(*,*) 'cont(',jcont,') = ',cont
2448:        write(*,*) 'theta = ',theta,'rhop = ',rhop
2449:        write(*,*) 'source = ',source,' sourceo = ',sourceo
2450:        write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2451:        write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2452:        write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2453:        write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2454:        write(*,*)
2455:       endif

2457:       return
2458:       end
2459:       double precision function  eexact(x,t)

2461:       implicit none

2463:       double precision x,t
2464:       double precision xot, head, tail, contact, ufan
2465:       double precision xpow, grat, urat
2466:       double precision uexact


2469:       logical debug

2471:       include 'ex74ftube.h'

2473:       debug = .false.


2476:       if (t .le. 0.0d+0) then
2477:         if (x .gt. 0.0d+0) then
2478:           eexact = e1
2479:         else
2480:           eexact = e4
2481:         endif
2482:       else

2484:        xot = x/t
2485:        head = -a4
2486:        tail = v3 - a3
2487:        contact = v2

2489:        if (xot .lt. head) then
2490:           eexact = e4
2491:        elseif (xot .gt. sspd) then
2492:           eexact = e1
2493:        elseif (xot .gt. contact) then
2494:           eexact = e2
2495:        elseif (xot .gt. tail) then
2496:           eexact = e3
2497:        else
2498:           ufan = uexact(x,t)
2499:           grat = (gamma - 1.0d+0) / 2.0d+0
2500:           xpow = 2.0d+0
2501:           urat = ufan / a4
2502:           eexact = e4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
2503:        endif

2505:       endif


2508:       if (debug) then
2509:         write(*,*)
2510:         write(*,*) 'eexact(',x,',',t,') = ',eexact
2511:         write(*,*)
2512:       endif

2514:       return
2515:       end
2516:       subroutine eigen(ht,uht)
2517: !23456789012345678901234567890123456789012345678901234567890123456789012
2518: !
2519: !          subroutine eigen
2520: !
2521: !  This subroutine computes the eigen values and eigen vectors
2522: !
2523: !23456789012345678901234567890123456789012345678901234567890123456789012


2526: !#######################################################################

2528:       implicit none

2530:       include 'ex74fcomd.h'

2532:       double precision ht, uht

2534:       double precision ut, at, lam1, lam2


2537: !#######################################################################

2539:       ut = uht / ht
2540:       at = sqrt( ht)

2542:       lam1 = ut - at
2543:       lam2 = ut + at

2545:       eigval(1) = lam1
2546:       eigval(2) = lam2

2548:       eigvec(1,1) = 1.0d+0
2549:       eigvec(2,1) = lam1
2550:       eigvec(1,2) = 1.0d+0
2551:       eigvec(2,2) = lam2

2553:       rinv(1,1) =  lam2 / (2.0d+0 * at)
2554:       rinv(2,1) = -lam1 / (2.0d+0 * at)
2555:       rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2556:       rinv(2,2) =  1.0d+0 / (2.0d+0 * at)


2559:       return
2560:       end
2561:       subroutine eigene(r,ru,e,l1,l2,l3)
2562: !23456789012345678901234567890123456789012345678901234567890123456789012
2563: !
2564: !          subroutine eigene
2565: !
2566: !  This subroutine computes the eigen values for the entropy fix
2567: !
2568: !23456789012345678901234567890123456789012345678901234567890123456789012


2571: !#######################################################################

2573:       implicit none

2575:       include 'ex74ftube.h'

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

2579:       double precision p,u,a

2581:       double precision eos
2582:       integer ierr

2584:       logical debug


2587: !#######################################################################

2589:       debug = .false.

2591:       if (debug) then
2592:          write(*,*)
2593:          write(*,*) 'gamma = ',gamma
2594:          write(*,*) 'r,ru,e = ',r,ru,e
2595:          write(*,*)
2596:       endif

2598:       p = eos(r,ru,e)
2599:       u = ru/r
2600:       if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2601:          write(*,*)
2602:          write(*,*) 'gamma = ',gamma
2603:          write(*,*) 'r = ',r
2604:          write(*,*) 'p = ',p
2605:          write(*,*)
2606:          call PetscFinalize(ierr)
2607:          stop
2608:       endif
2609:       a = sqrt((gamma * p)/r)

2611:       if (debug) then
2612:          write(*,*)
2613:          write(*,*) 'p,u,a = ',p,u,a
2614:          write(*,*)
2615:       endif

2617:       l1 = u - a
2618:       l2 = u
2619:       l3 = u + a

2621:       return
2622:       end
2623:       double precision function energy(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,jerg,xold)
2624: !
2625: !  This function computes the residual
2626: !  for the 1-D energy equation
2627: !
2628: !
2629:       implicit none

2631:       include 'ex74fcomd.h'
2632:       include 'ex74ftube.h'
2633: !
2634: !     input variables
2635: !
2636:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
2637:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
2638:       double precision ergee,  erge,  ergp,  ergw,  ergww
2639:       double precision xold(mx*neq)
2640: !
2641:       integer jerg
2642: !
2643: !     local variables
2644: !
2645:       double precision theta1
2646:       integer je
2647: !
2648: !  new
2649: !
2650:       double precision velfw, velfe
2651:       double precision vele,velp,velw
2652:       double precision fluxe, fluxw
2653:       double precision uepe, uepw
2654:       double precision ue, up, uw
2655:       double precision see, sep, sew
2656:       double precision seef, sewf
2657:       double precision upe, upw
2658:       double precision presse, pressw
2659:       double precision source
2660:       double precision te, tp, tw
2661:       double precision tef, twf, ref, rwf
2662:       double precision kef, kwf
2663:       double precision hflxe, hflxw
2664: !
2665: ! old
2666: !
2667:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
2668:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2669:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
2670:       double precision velfow, velfoe
2671:       double precision veloe,velop,velow
2672:       double precision fluxoe, fluxow
2673:       double precision uepoe, uepow
2674:       double precision uoe, uop, uow
2675:       double precision seoe, seop, seow
2676:       double precision seoef, seowf
2677:       double precision upoe, upow
2678:       double precision pressoe, pressow
2679:       double precision sourceo
2680:       double precision toe, top, tow
2681:       double precision toef, towf, roef, rowf
2682:       double precision koef, kowf
2683:       double precision hflxoe, hflxow
2684: !
2685: ! functions
2686: !
2687:       double precision godunov2, eos
2688:       double precision upwind, fluxlim

2690: !
2691: !
2692: ! ******************************************************************
2693: !
2694: !
2695: !
2696:       je = (neq*jerg)

2698: !########################
2699: !
2700: !      NEW
2701: !
2702: !########################

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

2706:       pressw  = eos(rhow, rhouw, ergw)
2707:       presse  = eos(rhoe, rhoue, erge)

2709:       uw = rhouw / rhow
2710:       up = rhoup / rhop
2711:       ue = rhoue / rhoe

2713:       upw = uw * pressw
2714:       upe = ue * presse

2716:       velfe = 0.5d+0 * (vele + velp)
2717:       velfw = 0.5d+0 * (velw + velp)

2719:       if (ihod .eq. 1) then

2721:         uepe = upwind(ergp,erge,velfe)
2722:         uepw = upwind(ergw,ergp,velfw)

2724:       elseif (ihod .eq. 2) then

2726:         uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2727:         uepw = fluxlim(ergww,ergw,ergp,erge,velfw)

2729:       endif

2731:       if (ihod .eq. 3) then
2732:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,3)
2733:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,3)
2734:       else
2735:         fluxe = (dt/dx) * ( uepe  + (0.5d+0*upe) )
2736:         fluxw = (dt/dx) * ( uepw  + (0.5d+0*upw) )
2737:       endif
2738: !
2739: !  radiation
2740: !
2741:       if (kappa0 .eq. 0.0d+0) then
2742:         source = 0.0d+0
2743:       else

2745:         see = (erge/rhoe) - (0.5d+0 * ue * ue)
2746:         sep = (ergp/rhop) - (0.5d+0 * up * up)
2747:         sew = (ergw/rhow) - (0.5d+0 * uw * uw)

2749:         seef = 0.5d+0 * (see + sep)
2750:         sewf = 0.5d+0 * (sew + sep)

2752:         te  = see / csubv
2753:         tp  = sep / csubv
2754:         tw  = sew / csubv

2756:         tef = seef / csubv
2757:         twf = sewf / csubv

2759:         ref = 0.5d+0 * (rhoe + rhop)
2760:         rwf = 0.5d+0 * (rhow + rhop)

2762:         kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2763:         kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)

2765:         if (wilson) then
2766:            kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2767:            kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2768:         endif

2770:         if ( debug .and. (kef .gt. 1.0d+10) ) then
2771:           write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2772:         endif
2773:         if ( debug .and. (kwf .gt. 1.0d+10) ) then
2774:           write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2775:         endif

2777:         hflxe = kef * (te - tp) / dx
2778:         hflxw = kwf * (tp - tw) / dx

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

2782:       endif

2784: !########################
2785: !
2786: !      OLD
2787: !
2788: !########################

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

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

2794:       pressow  = eos(rhoow, rhouow, ergow)
2795:       pressoe  = eos(rhooe, rhouoe, ergoe)

2797:       uow = rhouow / rhoow
2798:       uop = rhouop / rhoop
2799:       uoe = rhouoe / rhooe

2801:       upow = uow * pressow
2802:       upoe = uoe * pressoe

2804:       velfoe = 0.5d+0 * (veloe + velop)
2805:       velfow = 0.5d+0 * (velow + velop)


2808:       if (ihod .eq. 1) then

2810:         uepoe = upwind(ergop,ergoe,velfoe)
2811:         uepow = upwind(ergow,ergop,velfow)

2813:       elseif (ihod .eq. 2) then

2815:         uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
2816:         uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)

2818:       endif

2820:       if (ihod .eq. 3) then
2821:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,3)
2822:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,3)
2823:       else
2824:         fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
2825:         fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
2826:       endif

2828: !
2829: !  old radiation
2830: !
2831:       if (kappa0 .eq. 0.0d+0) then
2832:         sourceo = 0.0d+0
2833:       else

2835:         seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
2836:         seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
2837:         seow = (ergow/rhoow) - (0.5d+0 * uow * uow)

2839:         seoef = 0.5d+0 * (seoe + seop)
2840:         seowf = 0.5d+0 * (seow + seop)

2842:         toe  = seoe / csubv
2843:         top  = seop / csubv
2844:         tow  = seow / csubv

2846:         toef = seoef / csubv
2847:         towf = seowf / csubv

2849:         roef = 0.5d+0 * (rhooe + rhoop)
2850:         rowf = 0.5d+0 * (rhoow + rhoop)

2852:         koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
2853:         kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)

2855:         if (wilson) then
2856:            koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
2857:            kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
2858:         endif

2860:         if ( debug .and. (koef .gt. 1.0d+10) ) then
2861:           write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
2862:         endif
2863:         if ( debug .and. (kowf .gt. 1.0d+10) ) then
2864:           write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
2865:         endif

2867:         hflxoe = koef * (toe - top) / dx
2868:         hflxow = kowf * (top - tow) / dx

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

2872:       endif


2875: !########################
2876: !
2877: !      FUNCTION
2878: !
2879: !########################

2881: !VAM
2882:       if (probnum .eq. 3) then
2883:         fluxe  = 0.0d+0
2884:         fluxw  = 0.0d+0
2885:         fluxoe = 0.0d+0
2886:         fluxow = 0.0d+0
2887:       endif
2888: !VAM

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

2893:       if (debug) then
2894:         write(*,*)
2895:         write(*,*) 'energy(',jerg,') = ',energy
2896:         write(*,*)
2897:         write(*,*) fluxe,fluxw
2898:         write(*,*) fluxoe,fluxow
2899:         write(*,*) source,sourceo
2900:         write(*,*)
2901:       endif

2903:       return
2904:       end
2905:       double precision function eos(r,ru,e)

2907:       implicit none

2909:       double precision r,ru,e

2911:       double precision se, u

2913:       integer ierr

2915:       logical debug

2917:       include "ex74ftube.h"

2919:       debug = .false.

2921:       if (debug) then
2922:         write(*,*)
2923:         write(*,*) 'in eos r,ru,e'
2924:         write(*,*) r,ru,e
2925:         write(*,*)
2926:       endif

2928:       u = ru/r

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

2933:       if (eos .lt. 0.0d+0) then
2934:          write(*,*)
2935:          write(*,*) 'eos = ',eos
2936:          write(*,*) 'gamma = ',gamma
2937:          write(*,*) 'r = ',r
2938:          write(*,*) 'se = ',se
2939:          write(*,*) 'e = ',e
2940:          write(*,*) 'u = ',u
2941:          write(*,*) 'ru = ',ru
2942:          call PetscFinalize(ierr)
2943:          write(*,*)
2944:          stop
2945:       endif

2947:       if (debug) then
2948:         write(*,*)
2949:         write(*,*) 'in eos u,se,eos'
2950:         write(*,*) u,se,eos
2951:         write(*,*)
2952:       endif


2955:       return
2956:       end
2957:       subroutine eval2

2959:       implicit none

2961:       double precision prat, grat, xnum, xdenom


2964:       logical debug

2966:       include 'ex74ftube.h'

2968:       debug = .false.

2970:       prat = p2/p1
2971:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

2973:       xnum = grat + prat
2974:       xdenom = 1.0d+0 + (prat * grat)
2975: 
2976:       e2 = e1 * prat * (xnum/xdenom)
2977: 


2980:       if (debug) then
2981:         write(*,*)
2982:         write(*,*) 'e1  = ',e1
2983:         write(*,*) 'e2  = ',e2
2984:       endif

2986:       return
2987:       end
2988:       subroutine exact0

2990:       implicit none

2992:       double precision tol, xn
2993:       double precision shockp, fprime

2995:       integer maxnewt, niter

2997:       logical found, debug

2999:       include 'ex74ftube.h'

3001:       debug = .false.

3003:       tol = 1.0d-10

3005:       maxnewt = 40
3006: 
3007:       a1 = sqrt(gamma*p1/r1)
3008:       a4 = sqrt(gamma*p4/r4)



3012:       found = .false.
3013:       niter = 0

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

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

3019:         niter = niter + 1

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

3023:         if (debug) then
3024:           write(*,*) niter,shockp(xn),xn
3025:         endif

3027:         if ( abs(shockp(xn)) .lt. tol ) then
3028:            found = .true.
3029:         endif

3031:         goto 10

3033:       endif

3035:       if (.not. found) then

3037:          write(*,*) 'newton failed'
3038:          write(*,*) xn,shockp(xn)
3039:          stop

3041:       endif

3043:       p2 = xn


3046:       if (debug) then
3047:         write(*,*)
3048:         write(*,*) 'p1  = ',p1
3049:         write(*,*) 'p2  = ',p2
3050:         write(*,*) 'p4  = ',p4
3051:         write(*,*)
3052:       endif

3054:       return
3055:       end
3056:       double precision function flux(r,ru,e,eqn)
3057: !23456789012345678901234567890123456789012345678901234567890123456789012
3058: !
3059: !          function flux
3060: !
3061: !  This function computes the flux at a face
3062: !
3063: !23456789012345678901234567890123456789012345678901234567890123456789012


3066: !#######################################################################

3068:       implicit none

3070:       include 'ex74fcomd.h'
3071:       include 'ex74ftube.h'

3073:       double precision r, ru, e

3075:       integer eqn

3077:       double precision p,u

3079:       double precision eos


3082: !#######################################################################

3084:       p = eos(r,ru,e)
3085:       u = ru/r

3087:       if (eqn .eq. 1) then
3088:          flux = ru
3089:       elseif (eqn .eq. 2) then
3090:          flux = (u * ru) + p
3091:       else
3092:          flux = u * (e + p)
3093:       endif

3095:       return
3096:       end
3097:       double precision function fluxlim(fww,fw,fe,fee,vp)
3098: !23456789012345678901234567890123456789012345678901234567890123456789012
3099: !
3100: !          function fluxlim
3101: !
3102: !  this function computes the flux limited quick face value
3103: !
3104: !23456789012345678901234567890123456789012345678901234567890123456789012


3107: !#######################################################################

3109:       implicit none

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

3113:       double precision fd, fc, fu

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

3117:       double precision med, quick
3118: 
3119:       logical limit

3121: !#######################################################################

3123:       limit = .true.

3125:       if (vp .gt. 0.0d+0) then
3126:         fd = fe
3127:         fc = fw
3128:         fu = fww
3129:       else
3130:         fd = fw
3131:         fc = fe
3132:         fu = fee
3133:       endif

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

3137:       if (limit) then

3139:         beta = 0.25d+0
3140:         flc = 4.0d+0

3142:         f1 = fc
3143:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3144:         f3 = fu + ( flc * (fc - fu) )
3145:         f4 = med(f1,f2,f3)
3146:         fluxlim = vp * med(f1,f4,fhod)

3148:       else

3150:         fluxlim = vp * fhod

3152:       endif

3154:       return
3155:       end
3156:       double precision function fluxlim2(fww,fw,fe,fee,vp)
3157: !23456789012345678901234567890123456789012345678901234567890123456789012
3158: !
3159: !          function fluxlim2
3160: !
3161: !  this function computes the flux limited quick face value
3162: !
3163: !23456789012345678901234567890123456789012345678901234567890123456789012


3166: !#######################################################################

3168:       implicit none

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

3172:       double precision fd, fc, fu

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

3176:       double precision med, quick
3177: 
3178:       logical limit, debug

3180: !#######################################################################

3182:       debug = .false.

3184:       if (debug) then
3185:         write(*,*)
3186:         write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3187:         write(*,*) fee,fe,fw,fww
3188:         write(*,*)
3189:       endif

3191:       limit = .true.

3193:       if (vp .gt. 0.0d+0) then
3194:         fd = fe
3195:         fc = fw
3196:         fu = fww
3197:       else
3198:         fd = fw
3199:         fc = fe
3200:         fu = fee
3201:       endif

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

3205:       if (limit) then

3207:         beta = 0.25d+0
3208:         flc = 4.0d+0

3210:         f1 = fc
3211:         f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3212:         f3 = fu + ( flc * (fc - fu) )
3213:         f4 = med(f1,f2,f3)
3214:         fluxlim2 =  med(f1,f4,fhod)

3216:       else

3218:         fluxlim2 = fhod

3220:       endif

3222:       return
3223:       end
3224:       double precision function fprime(x)

3226:       implicit none

3228:       double precision  x, eps
3229:       double precision  shockp

3231:       eps = 1.0d-8

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

3235:       return
3236:       end
3237:       double precision function godunov2(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3238: !23456789012345678901234567890123456789012345678901234567890123456789012
3239: !
3240: !          function godunov2
3241: !
3242: !  this function computes the roe/godunov2 face value
3243: !
3244: !23456789012345678901234567890123456789012345678901234567890123456789012


3247: !#######################################################################

3249:       implicit none

3251:       include 'ex74fcomd.h'
3252:       include 'ex74ftube.h'

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

3256:       integer eqn

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

3260:       double precision hlle


3263: !#######################################################################

3265:       if (gorder .eq. 1) then
3266:         rrg  = rr
3267:         rlg  = rl
3268:         rurg = rur
3269:         rulg = rul
3270:         erg  = er
3271:         elg  = el
3272:       else
3273:         call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,rrg, rlg,rurg, rulg, erg, elg)
3274:       endif

3276: !VAM  if (ientro .eq. 0) then
3277: !VAM     godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3278: !VAM  elseif(ientro .eq. 1) then
3279: !VAM     godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3280: !VAM  else
3281:          godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3282: !VAM  endif


3285:       return
3286:       end
3287:       double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3288: !23456789012345678901234567890123456789012345678901234567890123456789012
3289: !
3290: !          function hlle
3291: !
3292: !  this function computes the roe/hlle face value
3293: !
3294: !23456789012345678901234567890123456789012345678901234567890123456789012


3297: !#######################################################################

3299:       implicit none

3301:       include 'ex74fcomd.h'
3302:       include 'ex74ftube.h'

3304:       double precision rrg,rlg,rurg,rulg,erg,elg
3305:       integer eqn

3307:       double precision laml1, laml2, laml3
3308:       double precision lamr1, lamr2, lamr3
3309:       double precision sl, sr


3312:       double precision flux

3314:       integer i, j, ispeed
3315: 

3317: !#######################################################################

3319:       ispeed = 1

3321:       do i = 1,neq
3322:         fr(i) = flux(rrg,rurg,erg,i)
3323:         fl(i) = flux(rlg,rulg,elg,i)
3324:       enddo

3326:       deltau(1) = rrg  - rlg
3327:       deltau(2) = rurg - rulg
3328:       deltau(3) = erg  - elg

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

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

3336: !VAM  if (ispeed .eq. 1) then
3337: !VAM    sl = min(laml1,lamt1)
3338: !VAM    sr = max(lamt2,lamr2)
3339: !VAM  else
3340:         sl = min(laml1,lamr1)
3341:         sr = max(laml3,lamr3)
3342: !VAM  endif


3345:       do i = 1,neq
3346:         froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )/(sr-sl)
3347:       enddo

3349:       hlle = froe(eqn)


3352:       return
3353:       end
3354:       double precision function med(x1,x2,x3)
3355: !23456789012345678901234567890123456789012345678901234567890123456789012
3356: !
3357: !          function med
3358: !
3359: !  this function computes the median of three numbers
3360: !
3361: !23456789012345678901234567890123456789012345678901234567890123456789012


3364: !#######################################################################

3366:       implicit none

3368:       double precision x1, x2, x3
3369:       double precision xhi, xlo

3371: !#######################################################################

3373:       xhi = max(x1,x2,x3)
3374:       xlo = min(x1,x2,x3)

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

3378:       return
3379:       end
3380:       double precision function mom(rhoee,  rhoe,  rhop,  rhow,  rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee,  erge,  ergp,  ergw,  ergww,jmom,xold)
3381: !
3382: !  This function computes the residual
3383: !  for the 1-D momentum equation
3384: !
3385: !
3386:       implicit none

3388:       include 'ex74fcomd.h'
3389:       include 'ex74ftube.h'
3390: !
3391: !     input variables
3392: !
3393:       double precision rhoee,  rhoe,  rhop,  rhow,  rhoww
3394:       double precision rhouee, rhoue, rhoup, rhouw, rhouww
3395:       double precision ergee,  erge,  ergp,  ergw,  ergww
3396:       double precision xold(mx*neq)
3397: !
3398:       integer jmom
3399: !
3400: !     local variables
3401: !
3402:       double precision theta1
3403:       integer jru
3404: !
3405: !  new
3406: !
3407:       double precision velfw, velfe
3408:       double precision vele,velp,velw
3409:       double precision fluxe, fluxw
3410:       double precision uurhoe, uurhow
3411:       double precision pressee, presse, pressp,pressw, pressww
3412:       double precision rupee, rupe, rupp, rupw, rupww
3413:       double precision uee, ue, up, uw, uww
3414:       double precision source
3415: !
3416: ! old
3417: !
3418:       double precision rhooee,  rhooe,  rhoop,  rhoow,  rhooww
3419:       double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3420:       double precision ergoee,  ergoe,  ergop,  ergow,  ergoww
3421:       double precision velfow, velfoe
3422:       double precision veloe,velop,velow
3423:       double precision fluxoe, fluxow
3424:       double precision uurhooe, uurhoow
3425:       double precision pressoee, pressoe, pressop, pressow, pressoww
3426:       double precision rupoee, rupoe, rupop, rupow, rupoww
3427:       double precision uoee, uoe, uop, uow, uoww
3428:       double precision sourceo

3430:       double precision eps
3431: !
3432: ! functions
3433: !
3434:       double precision godunov2, eos
3435:       double precision upwind, fluxlim
3436: !
3437: !
3438: ! ******************************************************************
3439: !
3440: !
3441:       eps = 1.0d-32
3442: !
3443:       jru = (neq*jmom) - 1

3445: !########################
3446: !
3447: !      NEW
3448: !
3449: !########################

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

3453:       presse  = eos(rhoe, rhoue, erge )
3454:       pressw  = eos(rhow, rhouw, ergw )

3456:       velfe = 0.5d+0 * (vele + velp)
3457:       velfw = 0.5d+0 * (velw + velp)

3459:       if (ihod .eq. 1) then

3461:         uurhoe = upwind(rhoup,rhoue,velfe)
3462:         uurhow = upwind(rhouw,rhoup,velfw)

3464:       elseif (ihod .eq. 2) then

3466:         uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3467:         uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)

3469:       endif

3471:       if (ihod .eq. 3) then
3472:         fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,2)
3473:         fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,2)
3474:       else
3475:         fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3476:         fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3477:       endif


3480:       source = 0.0d+0

3482: !########################
3483: !
3484: !      OLD
3485: !
3486: !########################

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

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

3492:       pressoe  = eos(rhooe, rhouoe, ergoe)
3493:       pressow  = eos(rhoow, rhouow, ergow)

3495:       velfoe = 0.5d+0 * (veloe + velop)
3496:       velfow = 0.5d+0 * (velow + velop)

3498:       if (ihod .eq. 1) then

3500:         uurhooe = upwind(rhouop,rhouoe,velfoe)
3501:         uurhoow = upwind(rhouow,rhouop,velfow)

3503:       elseif (ihod .eq. 2) then

3505:         uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3506:         uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)

3508:       endif

3510:       if (ihod .eq. 3) then
3511:         fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee,ergow, ergop, ergoe, ergoee,2)
3512:         fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,2)
3513:       else
3514:         fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3515:         fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3516:       endif

3518:       sourceo = 0.0d+0


3521: !########################
3522: !
3523: !      FUNCTION
3524: !
3525: !########################

3527:       theta1 = 1.0d+0 - theta
3528:       mom =  (rhoup - xold(jru))  + (  theta  * ( (fluxe  - fluxw ) - source  )  )  + (  theta1 * ( (fluxoe - fluxow) - sourceo )  )
3529: !VAM
3530:       if (probnum .eq. 3) then
3531:         mom = 0.0d+0
3532:       endif
3533: !VAM
3534:       if (debug) then
3535:         write(*,*)
3536:         write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3537:         write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3538:         write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3539:         write(*,*) 'presse = ',presse,'pressw = ',pressw
3540:         write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3541:         write(*,*)
3542:       endif

3544:       return
3545:       end
3546:       double precision function quick(fd, fc, fu)
3547: !23456789012345678901234567890123456789012345678901234567890123456789012
3548: !
3549: !          function quick
3550: !
3551: !  this function computes the quick face value
3552: !
3553: !23456789012345678901234567890123456789012345678901234567890123456789012


3556: !#######################################################################

3558:       implicit none

3560:       double precision fd, fc, fu

3562: !#######################################################################

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

3566:       return
3567:       end
3568:       double precision function  rexact(x,t)

3570:       implicit none

3572:       double precision x,t
3573:       double precision xot, head, tail, contact, ufan
3574:       double precision xpow, grat, urat
3575:       double precision uexact


3578:       logical debug

3580:       include 'ex74ftube.h'

3582:       debug = .false.


3585:       if (t .le. 0.0d+0) then
3586:         if (x .gt. 0.0d+0) then
3587:           rexact = r1
3588:         else
3589:           rexact = r4
3590:         endif
3591:       else

3593:        xot = x/t
3594:        head = -a4
3595:        tail = v3 - a3
3596:        contact = v2

3598:        if (xot .lt. head) then
3599:           rexact = r4
3600:        elseif (xot .gt. sspd) then
3601:           rexact = r1
3602:        elseif (xot .gt. contact) then
3603:           rexact = r2
3604:        elseif (xot .gt. tail) then
3605:           rexact = r3
3606:        else
3607:           ufan = uexact(x,t)
3608:           grat = (gamma - 1.0d+0) / 2.0d+0
3609:           xpow = 1.0d+0 / grat
3610:           urat = ufan / a4
3611:           rexact = r4 * (  ( 1.0d+0 - (grat * urat) ) ** xpow  )
3612:        endif

3614:       endif


3617:       if (debug) then
3618:         write(*,*)
3619:         write(*,*) 'rexact(',x,',',t,') = ',rexact
3620:         write(*,*)
3621:       endif

3623:       return
3624:       end
3625:       subroutine roestat(uhl, uhr, hl,hr,ht,uht)
3626: !23456789012345678901234567890123456789012345678901234567890123456789012
3627: !
3628: !          subroutine roestat
3629: !
3630: !  This subroutine computes the roe state at a face
3631: !
3632: !23456789012345678901234567890123456789012345678901234567890123456789012


3635: !#######################################################################

3637:       implicit none

3639:       include 'ex74fcomd.h'

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

3643:       double precision ul, ur, shl, shr, xnum, xdenom
3644: 


3647: !#######################################################################

3649:       ul = uhl / hl
3650:       ur = uhr / hr

3652:       shl = sqrt(hl)
3653:       shr = sqrt(hr)

3655:       xnum = (shl * ul) + (shr * ur)
3656:       xdenom = shl + shr

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

3661:       return
3662:       end
3663:       subroutine rval2

3665:       implicit none

3667:       double precision prat, grat, xnum, xdenom


3670:       logical debug

3672:       include 'ex74ftube.h'

3674:       debug = .false.

3676:       prat = p2/p1
3677:       grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)

3679:       xnum = 1.0d+0 + (grat * prat)
3680:       xdenom = grat + prat
3681: 
3682:       r2 = r1 * (xnum/xdenom)
3683: 


3686:       if (debug) then
3687:         write(*,*)
3688:         write(*,*) 'r1  = ',r1
3689:         write(*,*) 'r2  = ',r2
3690:       endif

3692:       return
3693:       end
3694:       subroutine  secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err, rrg, rlg,rurg, rulg, erg, elg)
3695: !23456789012345678901234567890123456789012345678901234567890123456789012
3696: !
3697: !          subroutine secondq
3698: !
3699: !  this subroutine computes the second order (based on quick) left
3700: !  and right states for the godunov solver.
3701: !
3702: !23456789012345678901234567890123456789012345678901234567890123456789012


3705: !#######################################################################

3707:       implicit none

3709:       include 'ex74fcomd.h'

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



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

3718:       double precision fluxlim2


3721: !#######################################################################

3723: !
3724: !  compute the velocities
3725: !
3726:       ull = rull/rll
3727:       ul  = rul /rl
3728:       ur  = rur /rr
3729:       urr = rurr/rrr

3731: !
3732: !  compute the left state first
3733: !
3734:       veld = 1.0d+0

3736:       rlg = fluxlim2(rll,rl,rr,rrr,veld)
3737:       ulg = fluxlim2(ull,ul,ur,urr,veld)
3738:       rulg = rlg * ulg
3739:       elg = fluxlim2(ell,el,er,err,veld)
3740: !
3741: !  now compute the right state
3742: !
3743:       veld = -1.0d+0

3745:       rrg = fluxlim2(rll,rl,rr,rrr,veld)
3746:       urg = fluxlim2(ull,ul,ur,urr,veld)
3747:       rurg = rrg * urg
3748:       erg = fluxlim2(ell,el,er,err,veld)



3752:       return
3753:       end
3754:       double precision function shockp(x)

3756:       implicit none

3758:       double precision x
3759:       double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
3760:       logical debug

3762:       include 'ex74ftube.h'

3764:       debug = .false.


3767:       if (debug) then
3768:          write(*,*)
3769:          write(*,*) 'gamma = ',gamma
3770:          write(*,*) 'a1 = ',a1
3771:          write(*,*) 'a4 = ',a4
3772:          write(*,*) 'p1 = ',p1
3773:          write(*,*) 'p2 = ',x
3774:          write(*,*)
3775:       endif

3777:       xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
3778:       xdenom = sqrt  (  2.0d+0 * gamma * ( (2.0d+0*gamma) + (gamma + 1.0d+0) * ((x/p1) - 1) )  )
3779:       xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)

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


3784:       if (debug) then
3785:          write(*,*)
3786:          write(*,*) 'xnum = ',xnum
3787:          write(*,*) 'gamma = ',gamma
3788:          write(*,*) 'a1 = ',a1
3789:          write(*,*) 'a4 = ',a4
3790:          write(*,*) 'p1 = ',p1
3791:          write(*,*) 'xdenom = ',xdenom
3792:          write(*,*) 'xpow = ',xpow
3793:          write(*,*) 'shockp = ',shockp
3794:          write(*,*) 'p2 = ',x
3795:          write(*,*)
3796:       endif

3798:       return
3799:       end
3800:       double precision function  uexact(x,t)

3802:       implicit none

3804:       double precision x,t
3805:       double precision xot, head, tail


3808:       logical debug

3810:       include 'ex74ftube.h'

3812:       debug = .false.

3814:       if (debug) then
3815:         write(*,*)
3816:         write(*,*) 't = ',t
3817:         write(*,*) 'x = ',x
3818:         write(*,*) 'a4 = ',a4
3819:         write(*,*) 'v3 = ',v3
3820:         write(*,*) 'a3 = ',a3
3821:         write(*,*)
3822:       endif

3824:       if (t .le. 0.0d+0) then
3825:         uexact = 0.0d+0
3826:       else

3828:        xot = x/t
3829:        head = -a4
3830:        tail = v3 - a3

3832:        if (xot .lt. head) then
3833:           uexact = 0.0d+0
3834:        elseif (xot .gt. sspd) then
3835:           uexact = 0.0d+0
3836:        elseif (xot .gt. tail) then
3837:           uexact = v2
3838:        else
3839:           uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
3840:        endif

3842:       endif


3845:       if (debug) then
3846:         write(*,*)
3847: !VAM    write(*,*) 'x = ',x,' t = ',t
3848:         write(*,*) 'uexact = ',uexact
3849:         write(*,*)
3850:       endif

3852:       return
3853:       end
3854:       double precision function upwind(fw, fe, vp)
3855: !23456789012345678901234567890123456789012345678901234567890123456789012
3856: !
3857: !          function upwind
3858: !
3859: !  this function computes the upwind face value
3860: !
3861: !23456789012345678901234567890123456789012345678901234567890123456789012


3864: !#######################################################################

3866:       implicit none

3868:       double precision fw, fe, vp

3870: !#######################################################################

3872:       if (vp .gt. 0.0) then
3873:          upwind = vp * fw
3874:       else
3875:          upwind = vp * fe
3876:       endif

3878:       return
3879:       end
3880:       subroutine uval2

3882:       implicit none

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


3887:       logical debug

3889:       include 'ex74ftube.h'

3891:       debug = .false.

3893:       prat = p2/p1
3894:       grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
3895:       grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
3896:       arat = a1/gamma

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

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

3902:       if (debug) then
3903:         write(*,*)
3904:         write(*,*) 'v2  = ',v2
3905:       endif

3907:       return
3908:       end
3909:       subroutine val3

3911:       implicit none

3913:       double precision prat, rpow, epow, p3t


3916:       logical debug

3918:       include 'ex74ftube.h'

3920:       debug = .false.


3923:       p3 = p2

3925:       prat = p3/p4

3927:       rpow = 1.0d+0 / gamma

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

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

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

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

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

3939:       if (debug) then
3940:         write(*,*)
3941:         write(*,*) 'a3 = ',a3
3942:         write(*,*) 'r3 = ',r3
3943:         write(*,*) 'e3 = ',e3
3944:         write(*,*) 'p3 = ',p3
3945:         write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
3946:         write(*,*)
3947:       endif

3949:       return
3950:       end
3951:       subroutine wval

3953:       implicit none

3955:       double precision prat, grat, xnum


3958:       logical debug

3960:       include 'ex74ftube.h'

3962:       debug = .false.

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

3967:       xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
3968: 
3969:       sspd = a1 * sqrt(xnum)
3970: 


3973:       if (debug) then
3974:         write(*,*)
3975:         write(*,*) 'sspd  = ',sspd
3976:       endif

3978:       return
3979:       end