Actual source code: ex33f.F
1: >
2: program radhyd
3: ! "$Id: ex4f.F,v 1.39 1999/03/10 19:29:25 Vince Mousseau $";
4: !
5: ! Description: This example solves a nonlinear system on 1 processor with SNES.
6: ! We solve the Euler equations in one dimension.
7: ! The command line options include:
8: ! -dt <dt>, where <dt> indicates time step
9: ! -mx <xg>, where <xg> = number of grid points in the x-direction
10: ! -nstep <nstep>, where <nstep> = number of time steps
11: ! -debug <ndb>, where <ndb> = 0) no debug 1) debug
12: ! -pcnew <npc>, where <npc> = 0) no preconditioning 1) rad preconditioning
13: ! -probnum <probnum>, where <probnum> = 1) cyclic Riesner 2) dam break
14: ! -ihod <ihod>, where <ihod> = 1) upwind 2) quick 3) godunov
15: ! -ientro <ientro>, where <ientro> = 0) basic 1) entropy fix 2) hlle
16: ! -theta <theta>, where <theta> = 0-1 0-explicit 1-implicit
17: ! -hnaught <hnaught>, where <hnaught> = height of left side
18: ! -hlo <hlo>, where <hlo> = hieght of right side
19: ! -ngraph <ngraph>, where <ngraph> = number of time steps between graphics
20: ! -damfac <damfac>, where <damfac> = fractional downward change in hight
21: ! -dampit <ndamp>, where <ndamp> = 1 turn damping on
22: ! -gorder <gorder>, where <gorder> = spatial oerder of godunov
23: !
24: !
25: !
26: --------------------------------------------------------------------------
27: !
28: ! Shock tube example
29: !
30: ! In this example the application context is a Fortran integer array:
31: ! ctx(1) = shell preconditioner pressure matrix contex
32: ! ctx(2) = semi implicit pressure matrix
33: ! ctx(4) = xold - old time values need for time advancement
34: ! ctx(5) = mx - number of control volumes
35: ! ctx(6) = N - total number of unknowns
36: !
37: !
38: --------------------------------------------------------------------------
40: implicit none
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Include files
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! The following include statements are generally used in SNES Fortran
47: ! programs:
48: ! petscsys.h - base PETSc routines
49: ! vec.h - vectors
50: ! mat.h - matrices
51: ! ksp.h - Krylov subspace methods
52: ! pc.h - preconditioners
53: ! snes.h - SNES interface
54: ! In addition, we need the following for use of PETSc drawing routines
55: ! draw.h - drawing routines
56: ! Other include statements may be needed if using additional PETSc
57: ! routines in a Fortran program, e.g.,
58: ! viewer.h - viewers
59: ! is.h - index sets
60: !
61: #include finclude/petscsys.h
62: #include finclude/petscvec.h
63: #include finclude/petscis.h
64: #include finclude/petscdraw.h
65: #include finclude/petscmat.h
66: #include finclude/petscksp.h
67: #include finclude/petscpc.h
68: #include finclude/petscsnes.h
69: #include finclude/petscviewer.h
71: #include "comd.h"
72: #include "tube.h"
74: !
75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: ! Variable declarations
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: !
79: ! Variables:
80: ! snes - nonlinear solver
81: ! x,r - solution, residual vectors
82: ! J - Jacobian matrix
83: ! its - iterations for convergence
84: ! matrix_free - flag - 1 indicates matrix-free version
85: ! dt - time step size
86: ! draw - drawing context
87: !
88: PetscFortranAddr ctx(6)
89: integer nx,ny
90: SNES snes
91: KSP ksp
92: PC pc
93: Vec x,r
94: PetscViewer view0,view1,view2,
95: 1 view3, view4
96: Mat Psemi
97: integer matrix_free, flg, N, ierr, ngraph
98: integer nstep, ndt, size, rank, i
99: integer its, lits, totits, totlits
100: integer ndb, npc, ndamp, nwilson, ndtcon
101: double precision plotim
102: ! logical pcnew
104: double precision krtol,katol,kdtol
105: double precision natol,nrtol,nstol
106: integer kmit,nmit,nmf
109: ! Note: Any user-defined Fortran routines (such as FormJacobian)
110: ! MUST be declared as external.
112: external FormFunction, FormInitialGuess,FormDt,
113: & PCRadSetUp, PCRadApply, FormGraph,FormDampit
114: double precision eos
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: ! Initialize program
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: open (unit=87,file='Dt.out',status='unknown')
122: c
123: c start PETSc
124: c
125: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
126: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
127: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
129: if (size .ne. 1) then
130: if (rank .eq. 0) then
131: write(6,*) 'This is a uniprocessor example only!'
132: endif
133: SETERRQ(1,' ',ierr)
134: endif
136: ! Initialize problem parameters
138: debug = .false.
139: dampit = .false.
140: wilson = .true.
141: dtcon = .true.
142: pcnew = .true.
143: dtmax = 1.0d+2
144: dtmin = 1.00d-12
145: dt = 1.0d-2
146: mx = 100
147: nstep = 50
148: matrix_free = 1
149: probnum = 1
150: gorder = 1
152: tfinal = 1.0d+0
153: tplot = 0.2d+0
154: dtgrow = 1.05d+0
155: tcscal = 0.5d+0
156: hcscal = 0.5d+0
158: ihod = 3
159: ientro = 1
160: theta = 1.00d+0
161: pi = 3.14159d+0
163: zero = 0.0
164: ngraph = 10
166: ndb = 0
167: npc = 1
169: damfac = 0.9d+0
171: gamma = 1.25d+0
172: csubv = 1.0d+0 / (gamma - 1.0d+0)
174: v1 = 0.0d+0
175: v4 = 0.0d+0
177: e1 = 1.0d+0
178: e4 = 1.0d+0
180: r1 = 1.0d+0
181: r4 = 2.0d+0
183: ru1 = r1 * v1
184: ru4 = r4 * v4
186: et1 = r1 * ( (0.5d+0 * v1 * v1) + e1 )
187: et4 = r4 * ( (0.5d+0 * v4 * v4) + e4 )
189: p1 = eos(r1,ru1,et1)
190: p4 = eos(r4,ru4,et4)
192: a1 = sqrt(gamma*p1/r1)
193: a4 = sqrt(gamma*p4/r4)
195: erg0 = 1.0d+2
196: kappa0 = 1.0d+0
197: kappaa = -2.0d+0
198: kappab = 13.0d+0 / 2.0d+0
199: CVAM kappab = 2.5d+0
201: c
202: c load the command line options
203: c
204: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
205: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
206: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nstep',nstep,flg,
207: & ierr)
208: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-debug',ndb,flg,
209: & ierr)
210: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-pcnew',npc,flg,
211: & ierr)
212: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ihod',ihod,flg,
213: & ierr)
214: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ientro',ientro,flg,
215: & ierr)
216: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',
217: & theta,flg,ierr)
218: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ngraph',ngraph,flg,
219: & ierr)
220: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
221: & '-damfac',damfac,flg,ierr)
222: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dampit',ndamp,flg,
223: & ierr)
224: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
225: & '-wilson',nwilson,flg,ierr)
226: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-gorder',gorder,flg,
227: & ierr)
228: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
229: & '-probnum',probnum,flg,ierr)
230: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
231: & '-kappa0',kappa0,flg,ierr)
232: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
233: & '-erg0',erg0,flg,ierr)
234: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dtcon',ndtcon,flg,
235: & ierr)
236: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
237: & '-tfinal',tfinal,flg,ierr)
238: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
239: & '-tplot',tplot,flg,ierr)
240: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
241: & '-dtgrow',dtgrow,flg,ierr)
242: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
243: & '-tcscal',tcscal,flg,ierr)
244: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
245: & '-hcscal',hcscal,flg,ierr)
246: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
247: & '-dtmax',dtmax,flg,ierr)
249: if (ndamp .eq. 1) then
250: dampit = .true.
251: endif
253: if (nwilson .eq. 0) then
254: wilson = .false.
255: endif
257: if (ndb .eq. 1) then
258: debug = .true.
259: endif
261: if (npc .eq. 0) then
262: pcnew = .false.
263: endif
265: if (ndtcon .eq. 0) then
266: dtcon = .false.
267: endif
269: CVAM if (dt .ge. dtmax .or. dt .le. dtmin) then
270: CVAM if (rank .eq. 0) write(6,*) 'DT is out of range'
271: CVAM SETERRA(1,0,' ')
272: CVAM endif
274: N = mx*neq
276: ctx(5) = mx
277: ctx(6) = N
279: if (debug) then
280: write(*,*) 'mx = ',mx
281: endif
285: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: ! Create nonlinear solver context
287: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
291: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
292: ! Create vector data structures; set function evaluation routine
293: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295: CALL MatCreate(PETSC_COMM_WORLD, ctx(2), ierr)
296: CALL MatSetSizes(ctx(2),PETSC_DECIDE,PETSC_DECIDE,mx,mx,ierr)
297: CALL MatSetFromOptions(ctx(2),ierr)
298: if (debug) then
299: call MatGetSize(ctx(2),nx,ny,ierr)
300: write(*,*) 'number of rows = ',nx,' number of col = ',ny
301: endif
302: c
303: c full size vectors
304: c
305: CALL VecCreate(PETSC_COMM_WORLD,x,ierr)
306: CALL VecSetSizes(x,PETSC_DECIDE,N,ierr)
307: CALL VecSetType(x,VECMPI,ierr)
308: call VecSetFromOptions(x,ierr)
309: call VecDuplicate(x,r,ierr)
310: call VecDuplicate(x,ctx(4),ierr)
311: c
312: c set grid
313: c
314: dx = 2.0d+0/dfloat(mx)
315: xl0 = -1.0d+0 -(0.5d+0 * dx)
317: if (debug) then
318: write(*,*) 'dx = ',dx
319: endif
320:
322: ! Set function evaluation routine and vector. Whenever the nonlinear
323: ! solver needs to evaluate the nonlinear function, it will call this
324: ! routine.
325: ! - Note that the final routine argument is the user-defined
326: ! context that provides application-specific data for the
327: ! function evaluation routine.
329: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
331: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332: ! Customize nonlinear solver; set runtime options
333: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type
336: <type>)
338: call SNESSetFromOptions(snes,ierr)
339: c
340: c set the line search function to damp the newton update.
341: c
342: ! if (dampit) then
343: ! call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
344: ! endif
345: c
346: c get the linear solver info from the nonlinear solver
347: c
349: call SNESGetKSP(snes,ksp,ierr)
350: call KSPGetPC(ksp,pc,ierr)
351: ! call KSPGetKSP(ksp,ksp1,ierr)
352: CVAM call KSPSetType(ksp,KSPPREONLY,ierr)
353: call KSPSetType(ksp,KSPGMRES,ierr)
355: call KSPGetTolerances(ksp,krtol,katol,kdtol,kmit,ierr)
356: call SNESGetTolerances(snes,natol,nrtol,nstol,nmit,nmf,ierr)
358: write(*,*)
359: write(*,*)
360: write(*,*) 'Linear solver'
361: write(*,*)
362: write(*,*) 'rtol = ',krtol
363: write(*,*) 'atol = ',katol
364: write(*,*) 'dtol = ',kdtol
365: write(*,*) 'maxits = ',kmit
366: write(*,*)
367: write(*,*)
368: write(*,*) 'Nonlinear solver'
369: write(*,*)
370: write(*,*) 'rtol = ',nrtol
371: write(*,*) 'atol = ',natol
372: write(*,*) 'stol = ',nstol
373: write(*,*) 'maxits = ',nmit
374: write(*,*) 'max func = ',nmf
375: write(*,*)
376: write(*,*)
378: c
379: c Build shell based preconditioner if flag set
380: c
381: if (pcnew) then
382: call PCSetType(pc,PCSHELL,ierr)
383: call PCShellSetContext(pc,ctx,ierr)
384: call PCShellSetSetUp(pc,PCRadSetUp,ierr)
385: call PCShellSetApply(pc,PCRadApply,ierr)
386: endif
388: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
390: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391: ! Evaluate initial guess; then solve nonlinear system.
392: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
393: c
394: c initial counters
395: c
396: time = 0.0d+0
397: plotim = 0.0d+0
398: totits = 0
399: totlits = 0
401: ! Note: The user should initialize the vector, x, with the initial guess
402: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
403: ! to employ an initial guess of zero, the user should explicitly set
404: ! this vector to zero by calling VecSet().
406: call FormInitialGuess(x,ierr)
407: c
408: c open a window to plot results
409: c
410: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
411: & 'density',0,0,300,300,view0,ierr)
412: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
413: & 'velocity',320,0,300,300,view1,ierr)
414: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
415: & 'total energy',640,0,300,300,view2,ierr)
416: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
417: & 'temperature',0,360,300,300,view3,ierr)
418: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
419: & 'pressure',320,360,300,300,view4,ierr)
420: c
421: c graph initial conditions
422: c
423: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
424: c
425: c copy x into xold
426: c
427: call VecCopy(x,ctx(4),ierr)
428: call FormDt(snes,x,ctx,ierr)
429: c
430: c################################
431: c
432: c TIME STEP LOOP BEGIN
433: c
434: c################################
435: c
436: ndt = 0
438: 10 if ( (ndt .le. nstep) .and. ((time + 1.0d-10) .lt. tfinal) ) then
440: if (debug) then
441: write(*,*)
442: write(*,*) 'start of time loop'
443: write(*,*)
444: write(*,*) 'ndt = ',ndt
445: write(*,*) 'nstep = ',nstep
446: write(*,*) 'time = ',time
447: write(*,*) 'tfinal = ',tfinal
448: write(*,*)
449: endif
451: ndt = ndt + 1
452: c
453: c increment time
454: c
455: time = time + dt
456: plotim = plotim + dt
457: c
458: c call the nonlinear solver
459: c
460: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
461: CALL SNESGetIterationNumber(snes,its,ierr)
462: c
463: c get the number of linear iterations used by the nonlinear solver
464: c
465: call SNESGetNumberLinearIterations(snes,lits,ierr)
467: if (debug) then
468: write(*,*) 'in radhyd ',ndt,'x'
469: call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
470: endif
471: c
472: c increment the counters
473: c
474: totits = totits + its
475: totlits = totlits + lits
476: c
477: c Compute new time step
478: c
479: call FormDt(snes,x,ctx,ierr)
480: c
481: c Draw contour plot of solution
482: c
483: if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
484:
485: plotim = plotim - tplot
488: if (rank .eq. 0) then
489: write(6,100) totits,totlits,ndt,dt,time
490: endif
491: 100 format('Newt = ',i7,' lin =',i7,' step =',i7,
492: & ' dt = ',e8.3,' time = ',e10.4)
493: c
494: c graph state conditions
495: c
496: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
498: endif
499: c
500: c copy x into xold
501: c
502: call VecCopy(x,ctx(4),ierr)
505: goto 10
507: endif
508: c
509: c################################
510: c
511: c TIME STEP LOOP END
512: c
513: c################################
514: c
516: c
517: c graph final conditions
518: c
519: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
522: write(*,*)
523: write(*,*)
524: write(*,*) 'total Newton iterations = ',totits
525: write(*,*) 'total linear iterations = ',totlits
526: write(*,*) 'Average Newton per time step = ',
527: & dble(totits)/dble(ndt)
528: write(*,*) 'Average Krylov per Newton = ',
529: & dble(totlits)/dble(totits)
530: write(*,*)
531: write(*,*)
533: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
534: ! Free work space. All PETSc objects should be destroyed when they
535: ! are no longer needed.
536: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
539: call MatDestroy(ctx(2),ierr)
540: call VecDestroy(x,ierr)
541: call VecDestroy(ctx(4),ierr)
542: call VecDestroy(r,ierr)
543: call SNESDestroy(snes,ierr)
544: call PetscViewerDestroy(view0,ierr)
545: call PetscViewerDestroy(view1,ierr)
546: call PetscViewerDestroy(view2,ierr)
547: call PetscViewerDestroy(view3,ierr)
548: call PetscViewerDestroy(view4,ierr)
550: call PCDestroy(ctx(1),ierr)
552: call PetscFinalize(ierr)
554: close(87)
556: stop
557: end
558: subroutine ApplicationDampit(x,deltx,w,ierr)
559: ! ---------------------------------------------------------------------
560: !
561: ! ApplicationDampit - Damps the newton update, called by
562: ! the higher level routine FormDampit().
563: !
564: ! Input Parameter:
565: ! x - current iterate
566: ! deltx - update
567: !
568: ! Output Parameters:
569: ! x - new iterate
570: ! ierr - error code
571: !
572: ! Notes:
573: ! This routine only damps the density. May want to add energy
574: ! in the future
575: !
577: implicit none
579: ! Common blocks:
580: #include "comd.h"
582: ! Input/output variables:
583: PetscScalar x(mx*neq),deltx(mx*neq),
584: 1 w(mx*neq)
585: integer ierr
587: ! Local variables:
588: double precision facmin, fac, newx, xmin, se, dse
589: double precision u,en,rn,run
590: integer i, jr, jru, je
592: 0
594: if (debug) then
595: write(*,*) 'begin damping'
596: do i = 1,mx*neq
597: write(*,*)'x(',i,') = ',x(i)
598: enddo
599: write(*,*)
600: do i = 1,mx*neq
601: write(*,*)'deltx(',i,') = ',deltx(i)
602: enddo
603: endif
605: facmin = 1.0d+0
606: c
607: c set the scale factor
608: c
609: do i=1,mx
610: c
611: c set pointers
612: c
613: jr = (neq*i) - 2
614: jru = (neq*i) - 1
615: je = (neq*i)
616: c
617: c make sure dencity stayes positive
618: c
619: newx = x(jr) - deltx(jr)
620: xmin = damfac * x(jr)
622: if (newx .lt. xmin) then
623: fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
624: if (fac .lt. facmin) then
625: if (debug) then
626: write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
627: endif
628: facmin = fac
629: endif
630: endif
631: c
632: c make sure Total energy stayes positive
633: c
634: newx = x(je) - deltx(je)
635: xmin = damfac * x(je)
637: if (newx .lt. xmin) then
638: fac = (1.0d+0 - damfac)*x(je)/deltx(je)
639: if (fac .lt. facmin) then
640: if (debug) then
641: write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
642: endif
643: facmin = fac
644: endif
645: endif
646: c
647: c make sure specific internal energy stayes positive
648: c
649:
650: u = x(jru)/x(jr)
651: se = (x(je)/x(jr)) - (0.5d+0 * u * u)
653: en = x(je) - deltx(je)
654: rn = x(jr) - deltx(jr)
655: run = x(jru) - deltx(jru)
657: dse = se - ( (en/rn) - (0.5d+0 * (run/rn) * (run/rn)) )
660: newx = se - dse
661: xmin = damfac * se
663: if (newx .lt. xmin) then
664: fac = (1.0d+0 - damfac) * se / dse
665: if (fac .lt. facmin) then
666: if (debug) then
667: write(*,*) 'se',i, damfac,facmin,fac,se,dse
668: endif
669: facmin = fac
670: endif
671: endif
673: enddo
674: c
675: c write out warning
676: c
677: if (facmin .lt. 1.0d+0) then
678: write(*,*) 'facmin = ',facmin, damfac,time
679: c
680: c scale the vector
681: c
682: do i=1,neq*mx
683: w(i) = x(i) - (facmin * deltx(i))
684: enddo
685: else
686: do i=1,neq*mx
687: w(i) = x(i) - deltx(i)
688: enddo
689: endif
691: if (debug) then
692: write(*,*) 'end damping'
693: do i = 1,mx*neq
694: write(*,*) 'w(',i,') = ',w(i)
695: enddo
696: endif
698: return
699: end
700: subroutine ApplicationDt(x,xold,ierr)
701: ! ---------------------------------------------------------------------
702: !
703: ! ApplicationDt - compute CFL numbers. Called by
704: ! the higher level routine FormDt().
705: !
706: ! Input Parameter:
707: ! x - local vector data
708: !
709: ! Output Parameters:
710: ! ierr - error code
711: !
712: ! Notes:
713: ! This routine uses standard Fortran-style computations over a 2-dim
714: array.
715: !
717: implicit none
719: ! Common blocks:
720: #include "comd.h"
721: #include "tube.h"
723: ! Input/output variables:
724: PetscScalar x(mx*neq), xold(mx*neq)
725: integer ierr
727: ! Local variables:
728: integer i, jr, jru, je
729: c
730: c new
731: c
732: double precision rhoee, rhoe, rhop, rhow, rhoww,
733: & rhouee, rhoue, rhoup, rhouw, rhouww,
734: & ergee, erge, ergp, ergw, ergww,
735: & vele, velp, velw
736: double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
737: double precision tcflg, hcflg, dtt, dth
738: double precision te, tp, tw
739: double precision ue, up, uw
740: double precision see, sep, sew
741: c
742: c old
743: c
744: double precision rhooee, rhooe, rhoop, rhoow, rhooww
745: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
746: double precision ergoee, ergoe, ergop, ergow, ergoww
747: double precision veloe, velop, velow
748: double precision uop, seop, top
749: double precision dtold, dttype
750: c
751: c functions
752: c
753: double precision eos
755: dtold = dt
757: 0
759: if (debug) then
760: write(*,*) 'in start dt'
761: do i = 1,mx*neq
762: write(*,*)'x(',i,') = ',x(i)
763: enddo
764: write(*,*) 'tfinal = ',tfinal
765: write(*,*) 'time = ',time
766: write(*,*) 'dt = ',dt
767: write(*,*) 'dtmax = ',dtmax
768: endif
770: sndp = -1.0d+20
771: vradn = 0.0d+0
772: vradd = 0.0d+0
774: c
775: c################################
776: c
777: c loop over all cells begin
778: c
779: c################################
780: c
781: do i=1,mx
782: c
783: c set pointers
784: c
785: jr = (neq*i) - 2
786: jru = (neq*i) - 1
787: je = (neq*i)
788: c
789: c
790: c set scalars
791: c
792: call Setpbc(i,x,
793: & rhoee, rhoe, rhop, rhow, rhoww,
794: & rhouee, rhoue, rhoup, rhouw, rhouww,
795: & ergee, erge, ergp, ergw, ergww,
796: & vele, velp, velw)
797: c
798: c compute temperatures
799: c
800: uw = rhouw / rhow
801: up = rhoup / rhop
802: ue = rhoue / rhoe
804: see = (erge/rhoe) - (0.5d+0 * ue * ue)
805: sep = (ergp/rhop) - (0.5d+0 * up * up)
806: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
808: te = see / csubv
809: tp = sep / csubv
810: tw = sew / csubv
811: c
812: c compute old temperature
813: c
815: call Setpbc(i,xold,
816: & rhooee, rhooe, rhoop, rhoow, rhooww,
817: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
818: & ergoee, ergoe, ergop, ergow, ergoww,
819: & veloe, velop, velow)
821: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
822: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
823: & ergoee, ergoe, ergop, ergow, ergoww,
824: & veloe, velop, velow, i)
826: uop = rhouop / rhoop
828: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
830: top = seop / csubv
831: c
832: c compute thermal cfl
833: c
834: vradn = vradn + (abs(tp - top)/dt)
835: vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
836: c
837: c compute hydro cfl
838: c
840: pressp = eos(rhop, rhoup, ergp)
841: sndp = max(sndp,sqrt( (gamma*pressp) / rhop ))
843: enddo
844: c
845: c################################
846: c
847: c loop over all cells end
848: c
849: c################################
850: c
852: vrad = vradn / vradd
854: tcfl = (vrad * dt) / dx
855: hcfl = (sndp * dt) / dx
857: dtt = max(dx/vrad,1.0d-7)
858: dtt = tcscal * dtt
860: dth = hcscal * dx / sndp
862: if (.not. dtcon) then
863: dt = min (dth,dtt,dt*dtgrow)
864: if (dt .lt. dtmin) then
865: dt = dtmin
866: endif
867: if (dt .gt. dtmax) then
868: dt = dtmax
869: endif
870: if ( (time + dt) .gt. tfinal) then
871: dt = tfinal - time
872: endif
874: if (dt .eq. dth) then
875: dttype = 1.0d+0
876: elseif (dt .eq. dtt) then
877: dttype = 2.0d+0
878: elseif (dt .eq. dtold*dtgrow) then
879: dttype = 3.0d+0
880: elseif (dt .eq. dtmax) then
881: dttype = 4.0d+0
882: elseif (dt .eq. dtmin) then
883: dttype = 5.0d+0
884: elseif (dt .eq. tfinal - time) then
885: dttype = 6.0
886: else
887: dttype = -1.0d+0
888: endif
890: endif
891:
892:
893: write(87,1000) time,dt,dth/hcscal,dtt/tcscal
894: write(88,1000) time,dttype
896: 1000 format(4(2x,e18.9))
898: if (debug) then
899: write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
900: write(*,*) 'dtt = ',dtt,' dth = ',dth
901: write(*,*) 'tfinal = ',tfinal
902: write(*,*) 'time = ',time
903: write(*,*) 'dt = ',dt
904: write(*,*) 'dtmax = ',dtmax
905: write(*,*)
906: endif
908: return
909: end
910: subroutine ApplicationExact(x,ierr)
911: ! ---------------------------------------------------------------------
912: !
913: ! ApplicationExact - Computes exact solution, called by
914: ! the higher level routine FormExact().
915: !
916: ! Input Parameter:
917: ! x - local vector data
918: !
919: ! Output Parameters:
920: ! x - initial conditions
921: ! ierr - error code
922: !
923: ! Notes:
924: ! This routine uses standard Fortran-style computations over a 1-dim
925: array.
926: !
928: implicit none
930: ! Common blocks:
932: #include "comd.h"
934: ! Input/output variables:
935: PetscScalar x(mx)
936: integer ierr
938: ! Local variables:
939: integer i
940: double precision xloc
941: PetscScalar rexact
944: ! Set parameters
946: 0
948: do i = 1,mx
950: xloc = xl0 + (dble(i) * dx)
951: x(i) = rexact(xloc,time)
953: enddo
955: return
956: end
957: subroutine ApplicationFunction(x,f,xold,ierr)
958: ! ---------------------------------------------------------------------
959: !
960: ! ApplicationFunction - Computes nonlinear function, called by
961: ! the higher level routine FormFunction().
962: !
963: ! Input Parameter:
964: ! x - local vector data
965: !
966: ! Output Parameters:
967: ! f - local vector data, f(x)
968: ! ierr - error code
969: !
970: ! Notes:
971: ! This routine uses standard Fortran-style computations over a 2-dim
972: array.
973: !
975: implicit none
977: ! Common blocks:
978: #include "comd.h"
980: ! Input/output variables:
981: PetscScalar x(mx*neq),f(mx*neq),
982: 1 xold(mx*neq)
983: integer ierr
985: ! Local variables:
986: integer i, jr, jru, je
987: double precision rhoee, rhoe, rhop, rhow, rhoww,
988: & rhouee, rhoue, rhoup, rhouw, rhouww,
989: & ergee, erge, ergp, ergw, ergww,
990: & vele, velp, velw
992: double precision cont, energy, mom
994: 0
996: if (debug) then
997: write(*,*) 'in function'
998: do i = 1,mx*neq
999: write(*,*)'x(',i,') = ',x(i)
1000: enddo
1001: endif
1002: c
1003: c################################
1004: c
1005: c loop over all cells begin
1006: c
1007: c################################
1008: c
1009: do i=1,mx
1010: c
1011: c set pointers
1012: c
1013: jr = (neq*i) - 2
1014: jru = (neq*i) - 1
1015: je = (neq*i)
1016: c
1017: c
1018: c set scalars
1019: c
1020: call Setpbc(i,x,
1021: & rhoee, rhoe, rhop, rhow, rhoww,
1022: & rhouee, rhoue, rhoup, rhouw, rhouww,
1023: & ergee, erge, ergp, ergw, ergww,
1024: & vele, velp, velw)
1025: c
1026: c compute functions
1027: c
1029: f(jr) = cont(rhoee, rhoe, rhop, rhow, rhoww,
1030: & rhouee, rhoue, rhoup, rhouw, rhouww,
1031: & ergee, erge, ergp, ergw, ergww,
1032: & i,xold)
1035: f(jru) = mom(rhoee, rhoe, rhop, rhow, rhoww,
1036: & rhouee, rhoue, rhoup, rhouw, rhouww,
1037: & ergee, erge, ergp, ergw, ergww,
1038: & i,xold)
1041: f(je) = energy(rhoee, rhoe, rhop, rhow, rhoww,
1042: & rhouee, rhoue, rhoup, rhouw, rhouww,
1043: & ergee, erge, ergp, ergw, ergww,
1044: & i,xold)
1046: if (debug) then
1047: write(*,*)
1048: write(*,*) i,jr,jru,je,'res,r,ru,e'
1049: write(*,*) f(jr),f(jru),f(je)
1050: write(*,*)
1051: endif
1053: enddo
1054: c
1055: c################################
1056: c
1057: c loop over all cells end
1058: c
1059: c################################
1060: c
1062: if (debug) then
1063: write(*,*) 'in function'
1064: do i = 1,mx*neq
1065: write(*,*) 'f(',i,') = ',f(i)
1066: enddo
1067: endif
1069: return
1070: end
1071: subroutine ApplicationInitialGuess(x,ierr)
1072: ! ---------------------------------------------------------------------
1073: !
1074: ! ApplicationInitialGuess - Computes initial approximation, called by
1075: ! the higher level routine FormInitialGuess().
1076: !
1077: ! Input Parameter:
1078: ! x - local vector data
1079: !
1080: ! Output Parameters:
1081: ! x - initial conditions
1082: ! ierr - error code
1083: !
1084: ! Notes:
1085: ! This routine uses standard Fortran-style computations over a 1-dim
1086: array.
1087: !
1089: implicit none
1091: ! Common blocks:
1093: #include "comd.h"
1094: #include "tube.h"
1096: ! Input/output variables:
1097: PetscScalar x(mx*neq)
1098: integer ierr
1100: ! Local variables:
1101: integer i, j, jr, jru, je
1102: double precision xloc, re, ee, ve
1103: double precision wid, efloor
1104: PetscScalar uexact, rexact, eexact
1107: CVAM efloor = max(1.0d-1,1.0d-3 * erg0)
1108: efloor = 1.0d-1
1109: CVAM wid = max(1.0d-2,dx)
1110: wid = 1.0d-2
1112: ! Set parameters
1114: 0
1116: do i = 1,mx
1118: jr = (neq*i) - 2
1119: jru = (neq*i) - 1
1120: je = (neq*i)
1122: xloc = xl0 + (dble(i) * dx)
1124: if (probnum .eq. 1) then
1125: re = rexact(xloc,zero)
1126: ve = uexact(xloc,zero)
1127: ee = eexact(xloc,zero)
1128: else
1129: re = 1.0d+0
1130: ve = 0.0d+0
1131: ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1132: endif
1134: x(jr) = re
1135: x(jru) = re * ve
1136: x(je) = re * ( (0.5d+0 * ve * ve) + ee )
1138: if (debug) then
1139: write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1140: 100 format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1141: endif
1143: enddo
1145: call exact0
1146: call eval2
1147: call rval2
1148: call wval
1149: call uval2
1150: v3 = v2
1151: call val3
1153: a1 = sqrt(gamma*p1/r1)
1154: a2 = sqrt(gamma*p2/r2)
1155: a3 = sqrt(gamma*p3/r3)
1156: a4 = sqrt(gamma*p4/r4)
1158: write(*,1000) r1,r2,r3,r4
1159: write(*,2000) p1,p2,p3,p4
1160: write(*,3000) e1,e2,e3,e4
1161: write(*,4000) a1,a2,a3,a4
1162: write(*,*)
1164: 1000 format ('rhos ',4(f12.6))
1165: 2000 format ('pressures ',4(f12.6))
1166: 3000 format ('energies ',4(f12.6))
1167: 4000 format ('sound ',4(f12.6))
1170: return
1171: end
1172: subroutine ApplicationXmgr(x,ivar,ierr)
1173: ! ---------------------------------------------------------------------
1174: !
1175: ! ApplicationXmgr - Sets the Xmgr output called from
1176: ! the higher level routine FormXmgr().
1177: !
1178: ! Input Parameter:
1179: ! x - local vector data
1180: !
1181: ! Output Parameters:
1182: ! x - initial conditions
1183: ! ierr - error code
1184: !
1185: ! Notes:
1186: ! This routine uses standard Fortran-style computations over a 1-dim
1187: array.
1188: !
1190: implicit none
1192: ! Common blocks:
1194: #include "comd.h"
1196: ! Input/output variables:
1197: PetscScalar x(mx)
1198: integer ivar,ierr
1200: ! Local variables:
1201: integer i
1202: double precision xloc, sum
1203: PetscScalar rexact
1204: integer iplotnum(5)
1205: save iplotnum
1206: character*8 grfile
1209: data iplotnum / -1,-1,-1,-1,-1 /
1213: ! Set parameters
1215: iplotnum(ivar) = iplotnum(ivar) + 1
1216: 0
1218: if (ivar .eq. 1) then
1219: write(grfile,4000) iplotnum(ivar)
1220: 4000 format('Xmgrr',i3.3)
1221: elseif (ivar .eq. 2) then
1222: write(grfile,5000) iplotnum(ivar)
1223: 5000 format('Xmgru',i3.3)
1224: elseif (ivar .eq. 3) then
1225: write(grfile,6000) iplotnum(ivar)
1226: 6000 format('Xmgre',i3.3)
1227: elseif (ivar .eq. 4) then
1228: write(grfile,7000) iplotnum(ivar)
1229: 7000 format('Xmgrt',i3.3)
1230: else
1231: write(grfile,8000) iplotnum(ivar)
1232: 8000 format('Xmgrp',i3.3)
1233: endif
1235: open (unit=44,file=grfile,status='unknown')
1237: do i = 1,mx
1239: xloc = xl0 + (dble(i) * dx)
1240: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1241: write(44,1000) xloc, x(i), rexact(xloc,time)
1242: else
1243: write(44,1000) xloc, x(i)
1244: endif
1246: enddo
1248: 1000 format(3(e18.12,2x))
1249: close(44)
1251: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1252: sum = 0.0d+0
1253: do i = 1,mx
1254: xloc = xl0 + (dble(i) * dx)
1255: sum = sum + (x(i) - rexact(xloc,time)) ** 2
1256: enddo
1257: sum = sqrt(sum)
1259: write(*,*)
1260: write(*,*) 'l2norm of the density error is',sum
1261: write(*,*)
1262: endif
1265: return
1266: end
1267: subroutine FormDampit(snes,ctx,x,f,g,y,w,
1268: & fnorm,ynorm,gnorm,flag,ierr)
1269: ! ---------------------------------------------------------------------
1270: !
1271: ! FormDampit - damps the Newton update
1272: !
1273: ! Input Parameters:
1274: ! snes - the SNES context
1275: ! x - current iterate
1276: ! f - residual evaluated at x
1277: ! y - search direction (containes new iterate on output)
1278: ! w - work vector
1279: ! fnorm - 2-norm of f
1280: !
1281: ! In this example the application context is a Fortran integer array:
1282: ! ctx(1) = shell preconditioner pressure matrix contex
1283: ! ctx(2) = semi implicit pressure matrix
1284: ! ctx(4) = xold - old time values need for time advancement
1285: ! ctx(5) = mx - number of control volumes
1286: ! ctx(6) = N - total number of unknowns
1287: !
1288: ! Output Parameter:
1289: ! g - residual evaluated at new iterate y
1290: ! y - new iterate (contains search direction on input
1291: ! gnorm - 2-norm of g
1292: ! ynorm - 2-norm of search length
1293: ! flag - set to 0 if the line search succeeds; -1 on failure
1294: !
1295: ! Notes:
1296: ! This routine serves as a wrapper for the lower-level routine
1297: ! "ApplicationDampit", where the actual computations are
1298: ! done using the standard Fortran style of treating the local
1299: ! vector data as a multidimensional array over the local mesh.
1300: ! This routine merely accesses the local vector data via
1301: ! VecGetArray() and VecRestoreArray().
1302: !
1303: implicit none
1305: #include finclude/petscsys.h
1306: #include finclude/petscvec.h
1307: #include finclude/petscsnes.h
1309: ! Input/output variables:
1310: SNES snes
1311: Vec x, f, g, y, w
1312: PetscFortranAddr ctx(*)
1313: PetscScalar fnorm, ynorm, gnorm
1314: integer ierr, flag
1316: ! Common blocks:
1318: #include "comd.h"
1320: ! Local variables:
1322: ! Declarations for use with local arrays:
1323: PetscScalar lx_v(0:1),ly_v(0:1),
1324: 1 lw_v(0:1)
1325: PetscOffset lx_i,ly_i,lw_i
1327: c
1328: c set ynorm
1329: c
1330: call VecNorm(y,NORM_2,ynorm,ierr)
1331: c
1332: c copy x to w
1333: c
1334: call VecCopy(x,w,ierr)
1335: c
1336: c get pointers to x, y, w
1337: c
1339: call VecGetArray(x,lx_v,lx_i,ierr)
1340: call VecGetArray(y,ly_v,ly_i,ierr)
1341: call VecGetArray(w,lw_v,lw_i,ierr)
1342: c
1343: c Compute Damping
1344: c
1345: call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1346: c
1347: c Restore vectors x, y, w
1348: c
1349: call VecRestoreArray(x,lx_v,lx_i,ierr)
1350: call VecRestoreArray(y,ly_v,ly_i,ierr)
1351: call VecRestoreArray(w,lw_v,lw_i,ierr)
1352: c
1353: c copy w to y
1354: c
1355: call VecCopy(w,y,ierr)
1356: c
1357: c compute new residual
1358: c
1359: call SNESComputeFunction(snes,y,g,ierr)
1360: call VecNorm(g,NORM_2,gnorm,ierr)
1361: flag = 0
1363: if (debug) then
1364: write(*,*) 'form damp ynorm = ',ynorm
1365: write(*,*) 'gnorm = ',gnorm
1366: endif
1368: return
1369: end
1370: subroutine FormDt(snes,x,ctx,ierr)
1371: ! ---------------------------------------------------------------------
1372: !
1373: ! FormDt - Compute CFL numbers
1374: !
1375: ! Input Parameters:
1376: ! snes - the SNES context
1377: ! x - input vector
1378: !
1379: ! In this example the application context is a Fortran integer array:
1380: ! ctx(1) = shell preconditioner pressure matrix contex
1381: ! ctx(2) = semi implicit pressure matrix
1382: ! ctx(4) = xold - old time values need for time advancement
1383: ! ctx(5) = mx - number of control volumes
1384: ! ctx(6) = N - total number of unknowns
1385: !
1386: !
1387: ! Notes:
1388: ! This routine serves as a wrapper for the lower-level routine
1389: ! "ApplicationDt", where the actual computations are
1390: ! done using the standard Fortran style of treating the local
1391: ! vector data as a multidimensional array over the local mesh.
1392: ! This routine merely accesses the local vector data via
1393: ! VecGetArray() and VecRestoreArray().
1394: !
1395: implicit none
1397: #include finclude/petscsys.h
1398: #include finclude/petscvec.h
1399: #include finclude/petscsnes.h
1401: ! Input/output variables:
1402: SNES snes
1403: Vec x
1404: PetscFortranAddr ctx(*)
1405: integer ierr
1407: ! Common blocks:
1409: #include "comd.h"
1411: ! Local variables:
1413: ! Declarations for use with local arrays:
1414: PetscScalar lx_v(0:1)
1415: PetscOffset lx_i
1416: PetscScalar lxold_v(0:1)
1417: PetscOffset lxold_i
1419: c
1420: c get pointers to x, xold
1421: c
1423: call VecGetArray(x,lx_v,lx_i,ierr)
1424: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1425: c
1426: c Compute function
1427: c
1428: call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1429: c
1430: c Restore vectors x, xold
1431: c
1432: call VecRestoreArray(x,lx_v,lx_i,ierr)
1433: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1435: return
1436: end
1437: subroutine FormExact(x,ierr)
1438: ! ---------------------------------------------------------------------
1439: !
1440: ! FormExact - Forms exact solution
1441: !
1442: ! Input Parameter:
1443: ! x - vector
1444: !
1445: ! Output Parameters:
1446: ! x - vector
1447: ! ierr - error code
1448: !
1449: ! Notes:
1450: ! This routine serves as a wrapper for the lower-level routine
1451: ! "ApplicationExact", where the actual computations are
1452: ! done using the standard Fortran style of treating the local
1453: ! vector data as a multidimensional array over the local mesh.
1454: ! This routine merely accesses the local vector data via
1455: ! VecGetArray() and VecRestoreArray().
1456: !
1457: implicit none
1459: #include finclude/petscsys.h
1460: #include finclude/petscvec.h
1461: #include finclude/petscmat.h
1462: #include finclude/petscsnes.h
1464: ! Input/output variables:
1465: Vec x
1466: integer ierr
1468: ! Declarations for use with local arrays:
1469: PetscScalar lx_v(0:1)
1470: PetscOffset lx_i
1472: 0
1474: c
1475: c get a pointer to x
1476: c
1477: call VecGetArray(x,lx_v,lx_i,ierr)
1478: c
1479: c Compute initial guess
1480: c
1481: call ApplicationExact(lx_v(lx_i),ierr)
1482: c
1483: c Restore vector x
1484: c
1485: call VecRestoreArray(x,lx_v,lx_i,ierr)
1487: return
1488: end
1489: subroutine FormFunction(snes,x,f,ctx,ierr)
1490: ! ---------------------------------------------------------------------
1491: !
1492: ! FormFunction - Evaluates nonlinear function, f(x).
1493: !
1494: ! Input Parameters:
1495: ! snes - the SNES context
1496: ! x - input vector
1497: !
1498: ! In this example the application context is a Fortran integer array:
1499: ! ctx(1) = shell preconditioner pressure matrix contex
1500: ! ctx(2) = semi implicit pressure matrix
1501: ! ctx(4) = xold - old time values need for time advancement
1502: ! ctx(5) = mx - number of control volumes
1503: ! ctx(6) = N - total number of unknowns
1504: !
1505: ! Output Parameter:
1506: ! f - vector with newly computed function
1507: !
1508: ! Notes:
1509: ! This routine serves as a wrapper for the lower-level routine
1510: ! "ApplicationFunction", where the actual computations are
1511: ! done using the standard Fortran style of treating the local
1512: ! vector data as a multidimensional array over the local mesh.
1513: ! This routine merely accesses the local vector data via
1514: ! VecGetArray() and VecRestoreArray().
1515: !
1516: implicit none
1518: #include finclude/petscsys.h
1519: #include finclude/petscvec.h
1520: #include finclude/petscsnes.h
1522: ! Input/output variables:
1523: SNES snes
1524: Vec x, f
1525: PetscFortranAddr ctx(*)
1526: integer ierr
1528: ! Common blocks:
1530: #include "comd.h"
1532: ! Local variables:
1534: ! Declarations for use with local arrays:
1535: PetscScalar lx_v(0:1), lf_v(0:1)
1536: PetscOffset lx_i, lf_i
1537: PetscScalar lxold_v(0:1)
1538: PetscOffset lxold_i
1540: c
1541: c get pointers to x, f, xold
1542: c
1544: call VecGetArray(x,lx_v,lx_i,ierr)
1545: call VecGetArray(f,lf_v,lf_i,ierr)
1546: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1547: c
1548: c Compute function
1549: c
1550: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),
1551: & lxold_v(lxold_i),ierr)
1552: c
1553: c Restore vectors x, f, xold
1554: c
1555: call VecRestoreArray(x,lx_v,lx_i,ierr)
1556: call VecRestoreArray(f,lf_v,lf_i,ierr)
1557: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1558: c
1559: c something to do with profiling
1560: c
1561: call PetscLogFlops(11.0d0*mx,ierr)
1563: return
1564: end
1565: subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1566: ! ---------------------------------------------------------------------
1567: !
1568: ! FormGraph - Forms Graphics output
1569: !
1570: ! Input Parameter:
1571: ! x - vector
1572: ! time - scalar
1573: !
1574: ! Output Parameters:
1575: ! ierr - error code
1576: !
1577: ! Notes:
1578: ! This routine serves as a wrapper for the lower-level routine
1579: ! "ApplicationXmgr", where the actual computations are
1580: ! done using the standard Fortran style of treating the local
1581: ! vector data as a multidimensional array over the local mesh.
1582: ! This routine merely accesses the local vector data via
1583: ! VecGetArray() and VecRestoreArray().
1584: !
1585: implicit none
1587: #include finclude/petscsys.h
1588: #include finclude/petscvec.h
1589: #include finclude/petscmat.h
1590: #include finclude/petscsnes.h
1592: #include "comd.h"
1593: #include "tube.h"
1595: ! Input/output variables:
1596: Vec x
1597: integer ierr
1599: ! Declarations for use with local arrays:
1600: IS rfrom,rto,
1601: 1 rufrom, ruto, efrom, eto
1602: Vec rval
1603: Vec uval
1604: Vec ruval
1605: Vec eval
1606: Vec seval
1607: Vec pval
1608: Vec kval
1609: Vec tval
1610: Vec steval
1611: VecScatter scatter
1612: PetscViewer view0,view1,
1613: 1 view2, view3, view4
1614: double precision minus1, l2err, gm1, csubvi
1617: csubvi = 1.0d+0 / csubv
1618: gm1 = gamma - 1.0d+0
1619: 0
1620: minus1 = -1.0d+0
1621: c
1622: c graphics vectors
1623: c
1624: CALL VecCreate(PETSC_COMM_WORLD,rval,ierr)
1625: CALL VecSetSizes(rval,PETSC_DECIDE,mx,ierr)
1626: CALL VecSetType(rval,VECMPI,ierr)
1627: call VecSetFromOptions(rval,ierr)
1628: call VecDuplicate(rval,uval,ierr)
1629: call VecDuplicate(rval,ruval,ierr)
1630: call VecDuplicate(rval,eval,ierr)
1631: call VecDuplicate(rval,seval,ierr)
1632: call VecDuplicate(rval,pval,ierr)
1633: call VecDuplicate(rval,kval,ierr)
1634: call VecDuplicate(rval,tval,ierr)
1635: call VecDuplicate(rval,steval,ierr)
1636: c
1637: c create index sets for vector scatters
1638: c
1639: call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1640: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, rto, ierr)
1641: call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1642: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, ruto, ierr)
1643: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1644: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, eto, ierr)
1646: c
1647: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1648: c
1649: c
1650: c load rval from x
1651: c
1652: call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1653: call VecScatterBegin(scatter,x,rval,INSERT_VALUES,
1654: & SCATTER_FORWARD,ierr)
1655: call VecScatterEnd(scatter,x,rval,INSERT_VALUES,
1656: & SCATTER_FORWARD,ierr)
1657: call VecScatterDestroy(scatter,ierr)
1658: c
1659: c plot rval vector
1660: c
1661: call VecView(rval,view0,ierr)
1662: c
1663: c make xmgr plot of rval
1664: c
1665: call FormXmgr(rval,1,ierr)
1666: c
1667: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1668: c
1669: c
1670: c load eval from x
1671: c
1672: call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1673: call VecScatterBegin(scatter,x,eval,INSERT_VALUES,
1674: & SCATTER_FORWARD,ierr)
1675: call VecScatterEnd(scatter,x,eval,INSERT_VALUES,
1676: & SCATTER_FORWARD,ierr)
1677: call VecScatterDestroy(scatter,ierr)
1678: c
1679: c plot eval vector
1680: c
1681: call VecView(eval,view2,ierr)
1682: c
1683: c make xmgr plot of eval
1684: c
1685: call FormXmgr(eval,3,ierr)
1686: c
1687: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1688: c
1690: c
1691: c load ruval from x
1692: c
1693: call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1694: call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,
1695: & SCATTER_FORWARD,ierr)
1696: call VecScatterEnd(scatter,x,ruval,INSERT_VALUES,
1697: & SCATTER_FORWARD,ierr)
1698: call VecScatterDestroy(scatter,ierr)
1699: c
1700: c create u = ru / r
1701: c
1702: call VecPointwiseDivide(ruval,rval,uval,ierr)
1703: c
1704: c plot uval vector
1705: c
1706: call VecView(uval,view1,ierr)
1707: c
1708: c make xmgr plot of uval
1709: c
1710: call FormXmgr(uval,2,ierr)
1712: c
1713: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1714: c
1716: call VecPointwiseMult(kval,uval,uval,ierr)
1717: call VecScale(kval,0.5d+0,ierr)
1719: call VecPointwiseDivide(steval,eval,rval,ierr)
1720: call VecWAXPY(seval,-1.0d+0,kval,steval,ierr)
1722: call VecCopy(seval,tval,ierr)
1723: call VecScale(tval,csubvi,ierr)
1725: c
1726: c plot tval vector
1727: c
1728: call VecView(tval,view3,ierr)
1729: c
1730: c make xmgr plot of tval
1731: c
1732: call FormXmgr(tval,4,ierr)
1734: c
1735: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736: c
1738: call VecPointwiseMult(rval,seval,pval,ierr)
1739: call VecScale(pval,gm1,ierr)
1740: c
1741: c plot pval vector
1742: c
1743: call VecView(pval,view4,ierr)
1744: c
1745: c make xmgr plot of pval
1746: c
1747: call FormXmgr(pval,5,ierr)
1748: c
1749: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1750: c
1756: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1757: ! Free work space. All PETSc objects should be destroyed when they
1758: ! are no longer needed.
1759: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1761: call VecDestroy(rval, ierr)
1762: call VecDestroy(uval, ierr)
1763: call VecDestroy(ruval,ierr)
1764: call VecDestroy(eval, ierr)
1765: call VecDestroy(seval, ierr)
1766: call VecDestroy(pval, ierr)
1767: call VecDestroy(kval, ierr)
1768: call VecDestroy(tval, ierr)
1769: call VecDestroy(steval, ierr)
1771: call ISDestroy(rfrom, ierr)
1772: call ISDestroy(rto, ierr)
1774: call ISDestroy(rufrom,ierr)
1775: call ISDestroy(ruto, ierr)
1777: call ISDestroy(efrom, ierr)
1778: call ISDestroy(eto, ierr)
1781: return
1782: end
1783: subroutine FormInitialGuess(x,ierr)
1784: ! ---------------------------------------------------------------------
1785: !
1786: ! FormInitialGuess - Forms initial approximation.
1787: !
1788: ! Input Parameter:
1789: ! x - vector
1790: !
1791: ! Output Parameters:
1792: ! x - vector
1793: ! ierr - error code
1794: !
1795: ! Notes:
1796: ! This routine serves as a wrapper for the lower-level routine
1797: ! "ApplicationInitialGuess", where the actual computations are
1798: ! done using the standard Fortran style of treating the local
1799: ! vector data as a multidimensional array over the local mesh.
1800: ! This routine merely accesses the local vector data via
1801: ! VecGetArray() and VecRestoreArray().
1802: !
1803: implicit none
1805: #include finclude/petscsys.h
1806: #include finclude/petscvec.h
1807: #include finclude/petscmat.h
1808: #include finclude/petscsnes.h
1810: ! Input/output variables:
1811: Vec x
1812: integer ierr
1814: ! Declarations for use with local arrays:
1815: PetscScalar lx_v(0:1)
1816: PetscOffset lx_i
1818: 0
1820: c
1821: c get a pointer to x
1822: c
1823: call VecGetArray(x,lx_v,lx_i,ierr)
1824: c
1825: c Compute initial guess
1826: c
1827: call ApplicationInitialGuess(lx_v(lx_i),ierr)
1828: c
1829: c Restore vector x
1830: c
1831: call VecRestoreArray(x,lx_v,lx_i,ierr)
1833: return
1834: end
1835: subroutine FormXmgr(x,ivar,ierr)
1836: ! ---------------------------------------------------------------------
1837: !
1838: ! FormXmgr - Forms Xmgr output
1839: !
1840: ! Input Parameter:
1841: ! x - vector
1842: !
1843: ! Output Parameters:
1844: ! x - vector
1845: ! ierr - error code
1846: !
1847: ! Notes:
1848: ! This routine serves as a wrapper for the lower-level routine
1849: ! "ApplicationXmgr", where the actual computations are
1850: ! done using the standard Fortran style of treating the local
1851: ! vector data as a multidimensional array over the local mesh.
1852: ! This routine merely accesses the local vector data via
1853: ! VecGetArray() and VecRestoreArray().
1854: !
1855: implicit none
1857: #include finclude/petscsys.h
1858: #include finclude/petscvec.h
1859: #include finclude/petscmat.h
1860: #include finclude/petscsnes.h
1862: ! Input/output variables:
1863: Vec x
1864: integer ivar,ierr
1866: ! Declarations for use with local arrays:
1867: PetscScalar lx_v(0:1)
1868: PetscOffset lx_i
1870: 0
1872: c
1873: c get a pointer to x
1874: c
1875: call VecGetArray(x,lx_v,lx_i,ierr)
1876: c
1877: c make the graph
1878: c
1879: call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1880: c
1881: c Restore vector x
1882: c
1883: call VecRestoreArray(x,lx_v,lx_i,ierr)
1885: return
1886: end
1887: subroutine PCRadApply(pc,x,y,ierr)
1888: ! -------------------------------------------------------------------
1889: !
1890: ! PCRadApply - This routine demonstrates the use of a
1891: ! user-provided preconditioner.
1892: !
1893: ! Input Parameters:
1894: ! pc - preconditioner object
1895: ! x - input vector
1896: ! In this example the shell preconditioner application context
1897: ! is a Fortran integer array:
1898: ! ctx(1) = shell preconditioner pressure matrix contex
1899: ! ctx(2) = semi implicit pressure matrix
1900: ! ctx(4) = xold - old time values need for time advancement
1901: ! ctx(5) = mx - number of control volumes
1902: ! ctx(6) = N - total number of unknowns
1903: !
1904: ! Output Parameters:
1905: ! y - preconditioned vector
1906: ! ierr - error code (nonzero if error has been detected)
1907: !
1908: ! Notes:
1909: ! This code implements the Jacobi preconditioner plus the
1910: ! SOR preconditioner
1911: !
1913: implicit none
1915: #include finclude/petscsys.h
1916: #include finclude/petscvec.h
1918: #include "comd.h"
1919: c
1920: c Input
1921: c
1922: PC pc
1923: Vec x, y
1924: integer ierr
1925: c
1926: c Local
1927: c
1928: PetscFortranAddr ctx(*)
1929: IS defrom, deto
1930: Vec de, rese
1931: VecScatter scatter
1932: PetscScalar lde_v(0:1),lrese_v(0:1)
1933: PetscOffset lde_i, lrese_i
1935: call PCShellGetContext(pc,ctx,ierr)
1937: c
1938: c Identity preconditioner
1939: c
1940: call VecCopy(x,y,ierr)
1941: c
1942: c if kappa0 not equal to zero then precondition the radiation diffusion
1943: c
1944: if (kappa0 .ne. 0.0d+0) then
1945:
1947: c
1948: c Create needed vectors
1949: c
1950: CALL VecCreate(PETSC_COMM_WORLD,de,ierr)
1951: CALL VecSetSizes(de,PETSC_DECIDE,mx,ierr)
1952: CALL VecSetType(de,VECMPI,ierr)
1953: call VecSetFromOptions(de,ierr)
1954: call VecDuplicate(de,rese,ierr)
1955: c
1956: c create index sets for scatters
1957: c
1958: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1959: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1960: c
1961: c load rese from x
1962: c
1963: call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1964: call VecScatterBegin(scatter,x,rese,INSERT_VALUES,
1965: & SCATTER_FORWARD,ierr)
1966: call VecScatterEnd(scatter,x,rese,INSERT_VALUES,
1967: & SCATTER_FORWARD,ierr)
1968: call VecScatterDestroy(scatter,ierr)
1969: c
1970: c apply preconditioner
1971: c
1972: call PCApply(ctx(1),rese,de,ierr)
1974: if (debug) then
1975: write(*,*) 'PCRadApply dh is'
1976: call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1977: endif
1978: c
1979: c load de into y
1980: c
1981: call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1982: call VecScatterBegin(scatter,de,y,INSERT_VALUES,
1983: & SCATTER_FORWARD,ierr)
1984: call VecScatterEnd(scatter,de,y,INSERT_VALUES,
1985: & SCATTER_FORWARD,ierr)
1986: call VecScatterDestroy(scatter,ierr)
1988: if (debug) then
1989: write(*,*) 'PCRadApply y is'
1990: call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1991: endif
1993: call VecDestroy(de,ierr)
1994: call VecDestroy(rese,ierr)
1996: call ISDestroy(defrom,ierr)
1997: call ISDestroy(deto,ierr)
1999: endif
2002: return
2003: end
2004: subroutine PCRadSetUp(pc,ierr)
2005: !
2006: ! PCRadSetUp - This routine sets up a user-defined
2007: ! preconditioner context.
2008: !
2009: ! Input Parameters:
2010: ! In this example the shell preconditioner application context
2011: ! is a Fortran integer array:
2012: ! ctx(1) = shell preconditioner pressure matrix contex
2013: ! ctx(2) = semi implicit pressure matrix
2014: ! ctx(4) = xold - old time values need for time advancement
2015: ! ctx(5) = mx - number of control volumes
2016: ! ctx(6) = N - total number of unknowns
2017: !
2018: ! Output Parameter:
2019: ! ierr - error code (nonzero if error has been detected)
2020: !
2021: ! Notes:
2022: ! In this example, we define the shell preconditioner to be Jacobi
2023: ! method. Thus, here we create a work vector for storing the reciprocal
2024: ! of the diagonal of the preconditioner matrix; this vector is then
2025: ! used within the routine PCRadApply().
2026: !
2028: implicit none
2030: #include finclude/petscsys.h
2031: #include finclude/petscvec.h
2032: #include finclude/petscmat.h
2034: #include "comd.h"
2035: c
2036: c Input
2037: c
2038: PC pc
2039: integer ierr
2040: c
2041: c Local
2042: c
2043: PetscFortranAddr ctx(*)
2044: Vec eold
2045:
2046: PetscScalar le_v(0:1)
2047: PetscOffset le_i
2049: call PCShellGetContext(pc,ctx,ierr)
2050:
2051: c
2052: c create vector
2053: c
2054: CALL VecCreate(PETSC_COMM_WORLD,eold,ierr)
2055: CALL VecSetSizes(eold,PETSC_DECIDE,mx,ierr)
2056: CALL VecSetType(eold,VECMPI,ierr)
2057: call VecSetFromOptions(eold,ierr)
2058: c
2059: c set up the matrix based on xold
2060: c
2061: call Setmat(ctx,ierr)
2062: c
2063: c set up the preconditioner
2064: c
2065: call PCDestroy(ctx(1),ierr)
2066: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
2067: CVAM call PCSetType(ctx(1),PCJACOBI,ierr)
2068: call PCSetType(ctx(1),PCLU,ierr)
2069: ! call PCSetVector(ctx(1),eold,ierr)
2070: call PCSetOperators(ctx(1),ctx(2),ctx(2),
2071: & DIFFERENT_NONZERO_PATTERN,ierr)
2072: call PCSetUp(ctx(1),ierr)
2074: call VecDestroy(eold,ierr)
2077: return
2078: end
2079: subroutine Setmat(ctx,ierr)
2081: implicit none
2083: #include finclude/petscsys.h
2084: #include finclude/petscvec.h
2085: #include finclude/petscmat.h
2086: ! Common blocks:
2087: #include "comd.h"
2088: #include "tube.h"
2090: ! Input/output variables:
2091: PetscFortranAddr ctx(*)
2092: integer ierr
2094: ! Local variables:
2095: PetscScalar lx_v(0:1)
2096: PetscOffset lx_i
2098: double precision xmult, himh, hiph, diag, upper, lower
2099: double precision hi, hip1, him1
2100: double precision
2101: & rhoee, rhoe, rhop, rhow, rhoww,
2102: & rhouee, rhoue, rhoup, rhouw, rhouww,
2103: & ergee, erge, ergp, ergw, ergww,
2104: & ue, up, uw
2105: double precision see, sep, sew, seef, sewf, tef, twf,
2106: & ref, rwf, kef, kwf, xmulte, xmultw
2107: c
2108: integer im, nx, ny
2109: c
2110: c get pointers to xold
2111: c
2112: call VecGetArray(ctx(4),lx_v,lx_i,ierr)
2113:
2115: c
2116: c############################
2117: c
2118: c loop over all cells begin
2119: c
2120: c############################
2121: c
2122: do im = 1,mx
2123: c
2124: c set scalars
2125: c
2126: call Setpbc(im,lx_v(lx_i),
2127: & rhoee, rhoe, rhop, rhow, rhoww,
2128: & rhouee, rhoue, rhoup, rhouw, rhouww,
2129: & ergee, erge, ergp, ergw, ergww,
2130: & ue, up, uw)
2131: c
2132: c set diffusion coefficients
2133: c
2134: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2135: sep = (ergp/rhop) - (0.5d+0 * up * up)
2136: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2138: seef = 0.5d+0 * (see + sep)
2139: sewf = 0.5d+0 * (sew + sep)
2141: tef = seef / csubv
2142: twf = sewf / csubv
2144: ref = 0.5d+0 * (rhoe + rhop)
2145: rwf = 0.5d+0 * (rhow + rhop)
2147: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2148: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2150: if (wilson) then
2151: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2152: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2153: endif
2154: c
2155: c set coefficients
2156: c
2157: xmult = dt / (dx * dx * csubv)
2159: xmulte = xmult * kef
2160: xmultw = xmult * kwf
2162: upper = -(xmulte / rhoe)
2163: lower = -(xmultw / rhow)
2165: diag = 1.0d+0 + ( (xmulte + xmultw) / rhop )
2167: c
2168: c load coefficients into the matrix
2169: c
2170: call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)
2172: if (im .eq. 1) then
2173: call MatSetValues(ctx(2),1,im-1,1,im ,upper,
2174: 1 INSERT_VALUES,ierr)
2175: elseif (im .eq. mx) then
2176: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2177: 1 INSERT_VALUES,ierr)
2178: else
2179: call MatSetValues(ctx(2),1,im-1,1,im ,upper,
2180: 1 INSERT_VALUES,ierr)
2181: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2182: 1 INSERT_VALUES,ierr)
2183: endif
2186: enddo
2187: c
2188: c############################
2189: c
2190: c loop over all cells end
2191: c
2192: c############################
2193: c
2194:
2195: c
2196: c final load of matrix
2197: c
2198: call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2199: call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2201: if (debug) then
2202: call MatGetSize(ctx(2),nx,ny,ierr)
2203: write(*,*) 'in setup nx = ',nx,' ny = ',ny
2204: call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2205: endif
2207: call VecRestoreArray (ctx(4),lx_v,lx_i,ierr)
2211: return
2212: end
2213: subroutine Setpbc(i,x,
2214: & rhoee, rhoe, rhop, rhow, rhoww,
2215: & rhouee, rhoue, rhoup, rhouw, rhouww,
2216: & ergee, erge, ergp, ergw, ergww,
2217: & vele, velp, velw)
2219: implicit none
2221: ! Common blocks:
2222: #include "comd.h"
2224: ! Input/output variables:
2225: PetscScalar x(mx*neq)
2226: integer i
2227: double precision rhoee, rhoe, rhop, rhow, rhoww
2228: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2229: double precision ergee, erge, ergp, ergw, ergww
2230: double precision vele, velp, velw
2232: ! Local variables:
2233: integer jr, jru, je
2235: c
2236: c set pointers
2237: c
2238: jr = (neq*i) - 2
2239: jru = (neq*i) - 1
2240: je = (neq*i)
2242: if (debug) then
2243: write(*,*)
2244: write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2245: write(*,*)
2246: endif
2247:
2248: if (i .eq. 1) then
2250: rhoee = x(jr+(2*neq))
2251: rhoe = x(jr+neq)
2252: rhop = x(jr)
2253: rhow = x(jr)
2254: rhoww = x(jr)
2256: rhouee = x(jru+(2*neq))
2257: rhoue = x(jru+neq)
2258: rhoup = x(jru)
2259: rhouw = x(jru)
2260: rhouww = x(jru)
2262: ergee = x(je+(2*neq))
2263: erge = x(je+neq)
2264: ergp = x(je)
2265: ergw = x(je)
2266: ergww = x(je)
2268: velw = 0.0d+0
2269: velp = rhoup/rhop
2270: vele = rhoue/rhoe
2272: elseif (i .eq. 2) then
2274: rhoee = x(jr+(2*neq))
2275: rhoe = x(jr+neq)
2276: rhop = x(jr)
2277: rhow = x(jr-neq)
2278: rhoww = x(jr-neq)
2280: rhouee = x(jru+(2*neq))
2281: rhoue = x(jru+neq)
2282: rhoup = x(jru)
2283: rhouw = x(jru-neq)
2284: rhouww = x(jru-neq)
2286: ergee = x(je+(2*neq))
2287: erge = x(je+neq)
2288: ergp = x(je)
2289: ergw = x(je-neq)
2290: ergww = x(je-neq)
2292: velw = rhouw/rhow
2293: velp = rhoup/rhop
2294: vele = rhoue/rhoe
2296: elseif (i .eq. mx-1) then
2298: rhoee = x(jr+neq)
2299: rhoe = x(jr+neq)
2300: rhop = x(jr)
2301: rhow = x(jr-neq)
2302: rhoww = x(jr-(2*neq))
2304: rhouee = x(jru+neq)
2305: rhoue = x(jru+neq)
2306: rhoup = x(jru)
2307: rhouw = x(jru-neq)
2308: rhouww = x(jru-(2*neq))
2310: ergee = x(je+neq)
2311: erge = x(je+neq)
2312: ergp = x(je)
2313: ergw = x(je-neq)
2314: ergww = x(je-(2*neq))
2316: velw = rhouw/rhow
2317: velp = rhoup/rhop
2318: vele = rhoue/rhoe
2320: elseif (i .eq. mx) then
2322: rhoee = x(jr)
2323: rhoe = x(jr)
2324: rhop = x(jr)
2325: rhow = x(jr-neq)
2326: rhoww = x(jr-(2*neq))
2328: rhouee = x(jru)
2329: rhoue = x(jru)
2330: rhoup = x(jru)
2331: rhouw = x(jru-neq)
2332: rhouww = x(jru-(2*neq))
2334: ergee = x(je)
2335: erge = x(je)
2336: ergp = x(je)
2337: ergw = x(je-neq)
2338: ergww = x(je-(2*neq))
2340: velw = rhouw/rhow
2341: velp = rhoup/rhop
2342: vele = 0.0d+0
2344: else
2346: rhoee = x(jr+(2*neq))
2347: rhoe = x(jr+neq)
2348: rhop = x(jr)
2349: rhow = x(jr-neq)
2350: rhoww = x(jr-(2*neq))
2352: rhouee = x(jru+(2*neq))
2353: rhoue = x(jru+neq)
2354: rhoup = x(jru)
2355: rhouw = x(jru-neq)
2356: rhouww = x(jru-(2*neq))
2358: ergee = x(je+(2*neq))
2359: erge = x(je+neq)
2360: ergp = x(je)
2361: ergw = x(je-neq)
2362: ergww = x(je-(2*neq))
2364: velw = rhouw/rhow
2365: velp = rhoup/rhop
2366: vele = rhoue/rhoe
2368: endif
2370: if (debug) then
2371: write(*,*)
2372: write(*,*) 'in Setpbc ',i,jr,jru,je
2373: write(*,*) 'mx = ',mx
2374: write(*,*)
2375: endif
2378: return
2379: end
2380: subroutine Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2381: & rhouee, rhoue, rhoup, rhouw, rhouww,
2382: & ergee, erge, ergp, ergw, ergww,
2383: & ue, up, uw, jbc)
2385: implicit none
2387: ! Common blocks:
2388: #include "comd.h"
2390: ! Input/output variables:
2391: integer jbc
2392: double precision rhoee, rhoe, rhop, rhow, rhoww
2393: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2394: double precision ergee, erge, ergp, ergw, ergww
2395: double precision ue, up, uw
2397: ! Local variables:
2399: if (jbc .eq. 1) then
2400: rhoww = rhop
2401: rhow = rhop
2402: rhouww = rhoup
2403: rhouw = rhoup
2404: ergww = ergp
2405: ergw = ergp
2406: uw = 0.0d+0
2407: elseif (jbc .eq. 2) then
2408: rhoww = rhow
2409: rhouww = rhouw
2410: ergww = ergw
2411: uw = rhouw / rhow
2412: else
2413: uw = rhouw / rhow
2414: endif
2416: if (jbc .eq. mx) then
2417: rhoee = rhop
2418: rhoe = rhop
2419: rhouee = rhoup
2420: rhoue = rhoup
2421: ergee = ergp
2422: erge = ergp
2423: ue = 0.0d+0
2424: elseif (jbc .eq. mx-1) then
2425: rhoee = rhoe
2426: rhouee = rhoue
2427: ergee = erge
2428: ue = rhoue / rhoe
2429: else
2430: ue = rhoue / rhoe
2431: endif
2433: up = rhoup / rhop
2435: if (debug) then
2436: write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2437: endif
2440: return
2441: end
2442: double precision function cont
2443: & (rhoee, rhoe, rhop, rhow, rhoww,
2444: & rhouee, rhoue, rhoup, rhouw, rhouww,
2445: & ergee, erge, ergp, ergw, ergww,
2446: & jcont,xold)
2447: c
2448: c This function computes the residual
2449: c for the 1-D continuity equation
2450: c
2451: c
2452: implicit none
2454: include 'comd.h'
2455: include 'tube.h'
2456: c
2457: c input variables
2458: c
2459: double precision rhoee, rhoe, rhop, rhow, rhoww
2460: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2461: double precision ergee, erge, ergp, ergw, ergww
2462: double precision xold(mx*neq)
2463: c
2464: integer jcont
2465: c
2466: c local variables
2467: c
2468: double precision theta1
2469: integer jr
2470: c
2471: c new
2472: c
2473: double precision velfw, velfe
2474: double precision vele,velp,velw
2475: double precision fluxe, fluxw
2476: double precision urhoe, urhow
2477: double precision source
2478: c
2479: c old
2480: c
2481: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2482: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2483: double precision ergoee, ergoe, ergop, ergow, ergoww
2484: double precision teoee, teoe, teop, teow, teoww,
2485: & uoe, uoee, uop, uow, uoww
2486: double precision velfow, velfoe
2487: double precision veloe,velop,velow
2488: double precision fluxoe, fluxow
2489: double precision urhooe, urhoow
2490: double precision sourceo
2491: c
2492: c functions
2493: c
2494: double precision godunov2
2495: double precision upwind, fluxlim
2496: c
2497: c
2498: c ******************************************************************
2499: c
2500: c
2501: c
2502: if (debug) then
2503: write(*,*)
2504: write(*,*) 'in cont',jcont,' ihod = ',ihod
2505: write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2506: write(*,*) 'rhop = ',rhop
2507: write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2508: write(*,*)
2509: endif
2511: jr = (neq*jcont) - 2
2513: c########################
2514: c
2515: c NEW
2516: c
2517: c########################
2519: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2520: & rhouee, rhoue, rhoup, rhouw, rhouww,
2521: & ergee, erge, ergp, ergw, ergww,
2522: & vele, velp, velw, jcont)
2524: velfe = 0.5d+0 * (vele + velp)
2525: velfw = 0.5d+0 * (velw + velp)
2527: if (ihod .eq. 1) then
2529: urhoe = upwind(rhop,rhoe,velfe)
2530: urhow = upwind(rhow,rhop,velfw)
2532: elseif (ihod .eq. 2) then
2534: urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2535: urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)
2537: endif
2539: if (ihod .eq. 3) then
2540: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2541: & rhouw,rhoup,rhoue,rhouee,
2542: & ergw, ergp, erge, ergee,1)
2543: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2544: & rhouww,rhouw,rhoup,rhoue,
2545: & ergww, ergw, ergp, erge,1)
2546: else
2547: fluxe = (dt/dx) * urhoe
2548: fluxw = (dt/dx) * urhow
2549: endif
2552: source = 0.0d+0
2554: c########################
2555: c
2556: c OLD
2557: c
2558: c########################
2560: call Setpbc(jcont,xold,
2561: & rhooee, rhooe, rhoop, rhoow, rhooww,
2562: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2563: & ergoee, ergoe, ergop, ergow, ergoww,
2564: & veloe, velop, velow)
2566: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
2567: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2568: & ergoee, ergoe, ergop, ergow, ergoww,
2569: & veloe, velop, velow, jcont)
2571: velfoe = 0.5d+0 * (veloe + velop)
2572: velfow = 0.5d+0 * (velow + velop)
2575: if (ihod .eq. 1) then
2577: urhooe = upwind(rhoop,rhooe,velfoe)
2578: urhoow = upwind(rhoow,rhoop,velfow)
2580: elseif (ihod .eq. 2) then
2582: urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2583: urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)
2585: endif
2587: if (ihod .eq. 3) then
2588: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
2589: & rhouow,rhouop,rhouoe,rhouoee,
2590: & ergow, ergop, ergoe, ergoee,1)
2591: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
2592: & rhouoww,rhouow,rhouop,rhouoe,
2593: & ergoww, ergow, ergop, ergoe,1)
2594: else
2595: fluxoe = (dt/dx) * urhooe
2596: fluxow = (dt/dx) * urhoow
2597: endif
2599: sourceo = 0.0d+0
2602: c########################
2603: c
2604: c FUNCTION
2605: c
2606: c########################
2608: theta1 = 1.0d+0 - theta
2609: cont = (rhop - xold(jr))
2610: & + ( theta * ( (fluxe - fluxw ) - source ) )
2611: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2612: CVAM
2613: if (probnum .eq. 3) then
2614: cont = 0.0d+0
2615: endif
2616: CVAM
2619: if (debug) then
2620: write(*,*)
2621: write(*,*) 'cont(',jcont,') = ',cont
2622: write(*,*) 'theta = ',theta,'rhop = ',rhop
2623: write(*,*) 'source = ',source,' sourceo = ',sourceo
2624: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2625: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2626: write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2627: write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2628: write(*,*)
2629: endif
2631: return
2632: end
2633: double precision function eexact(x,t)
2635: implicit none
2637: double precision x,t
2638: double precision xot, head, tail, contact, ufan
2639: double precision xpow, grat, urat
2640: double precision uexact
2643: logical debug
2645: include 'tube.h'
2647: debug = .false.
2650: if (t .le. 0.0d+0) then
2651: if (x .gt. 0.0d+0) then
2652: eexact = e1
2653: else
2654: eexact = e4
2655: endif
2656: else
2658: xot = x/t
2659: head = -a4
2660: tail = v3 - a3
2661: contact = v2
2663: if (xot .lt. head) then
2664: eexact = e4
2665: elseif (xot .gt. sspd) then
2666: eexact = e1
2667: elseif (xot .gt. contact) then
2668: eexact = e2
2669: elseif (xot .gt. tail) then
2670: eexact = e3
2671: else
2672: ufan = uexact(x,t)
2673: grat = (gamma - 1.0d+0) / 2.0d+0
2674: xpow = 2.0d+0
2675: urat = ufan / a4
2676: eexact = e4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
2677: endif
2679: endif
2682: if (debug) then
2683: write(*,*)
2684: write(*,*) 'eexact(',x,',',t,') = ',eexact
2685: write(*,*)
2686: endif
2688: return
2689: end
2690: subroutine eigen(ht,uht)
2691: c23456789012345678901234567890123456789012345678901234567890123456789012
2692: c
2693: c subroutine eigen
2694: c
2695: c This subroutine computes the eigen values and eigen vectors
2696: c
2697: c23456789012345678901234567890123456789012345678901234567890123456789012
2700: c#######################################################################
2702: implicit none
2704: include 'comd.h'
2706: double precision ht, uht
2708: double precision ut, at, lam1, lam2
2711: c#######################################################################
2713: ut = uht / ht
2714: at = sqrt( ht)
2716: lam1 = ut - at
2717: lam2 = ut + at
2719: eigval(1) = lam1
2720: eigval(2) = lam2
2722: eigvec(1,1) = 1.0d+0
2723: eigvec(2,1) = lam1
2724: eigvec(1,2) = 1.0d+0
2725: eigvec(2,2) = lam2
2727: rinv(1,1) = lam2 / (2.0d+0 * at)
2728: rinv(2,1) = -lam1 / (2.0d+0 * at)
2729: rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2730: rinv(2,2) = 1.0d+0 / (2.0d+0 * at)
2733: return
2734: end
2735: subroutine eigene(r,ru,e,l1,l2,l3)
2736: c23456789012345678901234567890123456789012345678901234567890123456789012
2737: c
2738: c subroutine eigene
2739: c
2740: c This subroutine computes the eigen values for the entropy fix
2741: c
2742: c23456789012345678901234567890123456789012345678901234567890123456789012
2745: c#######################################################################
2747: implicit none
2749: include 'tube.h'
2751: double precision r,ru,e,l1,l2,l3
2753: double precision p,u,a
2755: double precision eos
2756: integer ierr
2758: logical debug
2761: c#######################################################################
2763: debug = .false.
2765: if (debug) then
2766: write(*,*)
2767: write(*,*) 'gamma = ',gamma
2768: write(*,*) 'r,ru,e = ',r,ru,e
2769: write(*,*)
2770: endif
2772: p = eos(r,ru,e)
2773: u = ru/r
2774: if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2775: write(*,*)
2776: write(*,*) 'gamma = ',gamma
2777: write(*,*) 'r = ',r
2778: write(*,*) 'p = ',p
2779: write(*,*)
2780: call PetscFinalize(ierr)
2781: stop
2782: endif
2783: a = sqrt((gamma * p)/r)
2785: if (debug) then
2786: write(*,*)
2787: write(*,*) 'p,u,a = ',p,u,a
2788: write(*,*)
2789: endif
2791: l1 = u - a
2792: l2 = u
2793: l3 = u + a
2795: return
2796: end
2797: double precision function energy
2798: & (rhoee, rhoe, rhop, rhow, rhoww,
2799: & rhouee, rhoue, rhoup, rhouw, rhouww,
2800: & ergee, erge, ergp, ergw, ergww,
2801: & jerg,xold)
2802: c
2803: c This function computes the residual
2804: c for the 1-D energy equation
2805: c
2806: c
2807: implicit none
2809: include 'comd.h'
2810: include 'tube.h'
2811: c
2812: c input variables
2813: c
2814: double precision rhoee, rhoe, rhop, rhow, rhoww
2815: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2816: double precision ergee, erge, ergp, ergw, ergww
2817: double precision xold(mx*neq)
2818: c
2819: integer jerg
2820: c
2821: c local variables
2822: c
2823: double precision theta1
2824: integer je
2825: c
2826: c new
2827: c
2828: double precision velfw, velfe
2829: double precision vele,velp,velw
2830: double precision fluxe, fluxw
2831: double precision uepe, uepw
2832: double precision ue, up, uw
2833: double precision see, sep, sew
2834: double precision seef, sewf
2835: double precision upe, upw
2836: double precision presse, pressw
2837: double precision source
2838: double precision te, tp, tw
2839: double precision tef, twf, ref, rwf
2840: double precision kef, kwf
2841: double precision hflxe, hflxw
2842: c
2843: c old
2844: c
2845: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2846: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2847: double precision ergoee, ergoe, ergop, ergow, ergoww
2848: double precision velfow, velfoe
2849: double precision veloe,velop,velow
2850: double precision fluxoe, fluxow
2851: double precision uepoe, uepow
2852: double precision uoe, uop, uow
2853: double precision seoe, seop, seow
2854: double precision seoef, seowf
2855: double precision upoe, upow
2856: double precision pressoe, pressow
2857: double precision sourceo
2858: double precision toe, top, tow
2859: double precision toef, towf, roef, rowf
2860: double precision koef, kowf
2861: double precision hflxoe, hflxow
2862: c
2863: c functions
2864: c
2865: double precision godunov2, eos
2866: double precision upwind, fluxlim
2868: c
2869: c
2870: c ******************************************************************
2871: c
2872: c
2873: c
2874: je = (neq*jerg)
2876: c########################
2877: c
2878: c NEW
2879: c
2880: c########################
2882: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2883: & rhouee, rhoue, rhoup, rhouw, rhouww,
2884: & ergee, erge, ergp, ergw, ergww,
2885: & vele, velp, velw, jerg)
2887: pressw = eos(rhow, rhouw, ergw)
2888: presse = eos(rhoe, rhoue, erge)
2890: uw = rhouw / rhow
2891: up = rhoup / rhop
2892: ue = rhoue / rhoe
2894: upw = uw * pressw
2895: upe = ue * presse
2897: velfe = 0.5d+0 * (vele + velp)
2898: velfw = 0.5d+0 * (velw + velp)
2900: if (ihod .eq. 1) then
2902: uepe = upwind(ergp,erge,velfe)
2903: uepw = upwind(ergw,ergp,velfw)
2905: elseif (ihod .eq. 2) then
2907: uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2908: uepw = fluxlim(ergww,ergw,ergp,erge,velfw)
2910: endif
2912: if (ihod .eq. 3) then
2913: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2914: & rhouw,rhoup,rhoue,rhouee,
2915: & ergw, ergp, erge, ergee,3)
2916: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2917: & rhouww,rhouw,rhoup,rhoue,
2918: & ergww, ergw, ergp, erge,3)
2919: else
2920: fluxe = (dt/dx) * ( uepe + (0.5d+0*upe) )
2921: fluxw = (dt/dx) * ( uepw + (0.5d+0*upw) )
2922: endif
2923: c
2924: c radiation
2925: c
2926: if (kappa0 .eq. 0.0d+0) then
2927: source = 0.0d+0
2928: else
2930: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2931: sep = (ergp/rhop) - (0.5d+0 * up * up)
2932: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2934: seef = 0.5d+0 * (see + sep)
2935: sewf = 0.5d+0 * (sew + sep)
2937: te = see / csubv
2938: tp = sep / csubv
2939: tw = sew / csubv
2941: tef = seef / csubv
2942: twf = sewf / csubv
2944: ref = 0.5d+0 * (rhoe + rhop)
2945: rwf = 0.5d+0 * (rhow + rhop)
2947: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2948: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2950: if (wilson) then
2951: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2952: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2953: endif
2955: if ( debug .and. (kef .gt. 1.0d+10) ) then
2956: write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2957: endif
2958: if ( debug .and. (kwf .gt. 1.0d+10) ) then
2959: write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2960: endif
2962: hflxe = kef * (te - tp) / dx
2963: hflxw = kwf * (tp - tw) / dx
2965: source = (dt/dx) * (hflxe - hflxw)
2967: endif
2969: c########################
2970: c
2971: c OLD
2972: c
2973: c########################
2975: call Setpbc(jerg,xold,
2976: & rhooee, rhooe, rhoop, rhoow, rhooww,
2977: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2978: & ergoee, ergoe, ergop, ergow, ergoww,
2979: & veloe, velop, velow)
2981: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
2982: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2983: & ergoee, ergoe, ergop, ergow, ergoww,
2984: & veloe, velop, velow, jerg)
2986: pressow = eos(rhoow, rhouow, ergow)
2987: pressoe = eos(rhooe, rhouoe, ergoe)
2989: uow = rhouow / rhoow
2990: uop = rhouop / rhoop
2991: uoe = rhouoe / rhooe
2993: upow = uow * pressow
2994: upoe = uoe * pressoe
2996: velfoe = 0.5d+0 * (veloe + velop)
2997: velfow = 0.5d+0 * (velow + velop)
3000: if (ihod .eq. 1) then
3002: uepoe = upwind(ergop,ergoe,velfoe)
3003: uepow = upwind(ergow,ergop,velfow)
3005: elseif (ihod .eq. 2) then
3007: uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
3008: uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)
3010: endif
3012: if (ihod .eq. 3) then
3013: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3014: & rhouow,rhouop,rhouoe,rhouoee,
3015: & ergow, ergop, ergoe, ergoee,3)
3016: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3017: & rhouoww,rhouow,rhouop,rhouoe,
3018: & ergoww, ergow, ergop, ergoe,3)
3019: else
3020: fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
3021: fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
3022: endif
3024: c
3025: c old radiation
3026: c
3027: if (kappa0 .eq. 0.0d+0) then
3028: sourceo = 0.0d+0
3029: else
3031: seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
3032: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
3033: seow = (ergow/rhoow) - (0.5d+0 * uow * uow)
3035: seoef = 0.5d+0 * (seoe + seop)
3036: seowf = 0.5d+0 * (seow + seop)
3038: toe = seoe / csubv
3039: top = seop / csubv
3040: tow = seow / csubv
3042: toef = seoef / csubv
3043: towf = seowf / csubv
3045: roef = 0.5d+0 * (rhooe + rhoop)
3046: rowf = 0.5d+0 * (rhoow + rhoop)
3048: koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
3049: kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)
3051: if (wilson) then
3052: koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
3053: kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
3054: endif
3056: if ( debug .and. (koef .gt. 1.0d+10) ) then
3057: write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
3058: endif
3059: if ( debug .and. (kowf .gt. 1.0d+10) ) then
3060: write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
3061: endif
3063: hflxoe = koef * (toe - top) / dx
3064: hflxow = kowf * (top - tow) / dx
3066: sourceo = (dt/dx) * (hflxoe - hflxow)
3068: endif
3071: c########################
3072: c
3073: c FUNCTION
3074: c
3075: c########################
3077: CVAM
3078: if (probnum .eq. 3) then
3079: fluxe = 0.0d+0
3080: fluxw = 0.0d+0
3081: fluxoe = 0.0d+0
3082: fluxow = 0.0d+0
3083: endif
3084: CVAM
3086: theta1 = 1.0d+0 - theta
3087: energy = (ergp - xold(je))
3088: & + ( theta * ( (fluxe - fluxw ) - source ) )
3089: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3091: if (debug) then
3092: write(*,*)
3093: write(*,*) 'energy(',jerg,') = ',energy
3094: write(*,*)
3095: write(*,*) fluxe,fluxw
3096: write(*,*) fluxoe,fluxow
3097: write(*,*) source,sourceo
3098: write(*,*)
3099: endif
3101: return
3102: end
3103: double precision function eos(r,ru,e)
3105: implicit none
3107: double precision r,ru,e
3109: double precision se, u
3111: integer ierr
3113: logical debug
3115: include "tube.h"
3117: debug = .false.
3119: if (debug) then
3120: write(*,*)
3121: write(*,*) 'in eos r,ru,e'
3122: write(*,*) r,ru,e
3123: write(*,*)
3124: endif
3126: u = ru/r
3128: se = (e/r) - (0.5d+0 * u * u)
3129: eos = (gamma - 1.0d+0) * r * se
3131: if (eos .lt. 0.0d+0) then
3132: write(*,*)
3133: write(*,*) 'eos = ',eos
3134: write(*,*) 'gamma = ',gamma
3135: write(*,*) 'r = ',r
3136: write(*,*) 'se = ',se
3137: write(*,*) 'e = ',e
3138: write(*,*) 'u = ',u
3139: write(*,*) 'ru = ',ru
3140: call PetscFinalize(ierr)
3141: write(*,*)
3142: stop
3143: endif
3145: if (debug) then
3146: write(*,*)
3147: write(*,*) 'in eos u,se,eos'
3148: write(*,*) u,se,eos
3149: write(*,*)
3150: endif
3153: return
3154: end
3155: subroutine eval2
3157: implicit none
3159: double precision prat, grat, xnum, xdenom
3162: logical debug
3164: include 'tube.h'
3166: debug = .false.
3168: prat = p2/p1
3169: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
3171: xnum = grat + prat
3172: xdenom = 1.0d+0 + (prat * grat)
3173:
3174: e2 = e1 * prat * (xnum/xdenom)
3175:
3178: if (debug) then
3179: write(*,*)
3180: write(*,*) 'e1 = ',e1
3181: write(*,*) 'e2 = ',e2
3182: endif
3184: return
3185: end
3186: subroutine exact0
3188: implicit none
3190: double precision tol, xn
3191: double precision shockp, fprime
3193: integer maxnewt, niter
3195: logical found, debug
3197: include 'tube.h'
3199: debug = .false.
3201: tol = 1.0d-10
3203: maxnewt = 40
3204:
3205: a1 = sqrt(gamma*p1/r1)
3206: a4 = sqrt(gamma*p4/r4)
3210: found = .false.
3211: niter = 0
3213: xn = 0.5d+0 * (p1 + p4)
3215: 10 if ( (.not. found) .and. (niter .le. maxnewt) ) then
3217: niter = niter + 1
3219: xn = xn - (shockp(xn) / fprime(xn))
3221: if (debug) then
3222: write(*,*) niter,shockp(xn),xn
3223: endif
3225: if ( abs(shockp(xn)) .lt. tol ) then
3226: found = .true.
3227: endif
3229: goto 10
3231: endif
3233: if (.not. found) then
3235: write(*,*) 'newton failed'
3236: write(*,*) xn,shockp(xn)
3237: stop
3239: endif
3241: p2 = xn
3244: if (debug) then
3245: write(*,*)
3246: write(*,*) 'p1 = ',p1
3247: write(*,*) 'p2 = ',p2
3248: write(*,*) 'p4 = ',p4
3249: write(*,*)
3250: endif
3252: return
3253: end
3254: double precision function flux(r,ru,e,eqn)
3255: c23456789012345678901234567890123456789012345678901234567890123456789012
3256: c
3257: c function flux
3258: c
3259: c This function computes the flux at a face
3260: c
3261: c23456789012345678901234567890123456789012345678901234567890123456789012
3264: c#######################################################################
3266: implicit none
3268: include 'comd.h'
3269: include 'tube.h'
3271: double precision r, ru, e
3273: integer eqn
3275: double precision p,u
3277: double precision eos
3280: c#######################################################################
3282: p = eos(r,ru,e)
3283: u = ru/r
3285: if (eqn .eq. 1) then
3286: flux = ru
3287: elseif (eqn .eq. 2) then
3288: flux = (u * ru) + p
3289: else
3290: flux = u * (e + p)
3291: endif
3293: return
3294: end
3295: double precision function fluxlim(fww,fw,fe,fee,vp)
3296: c23456789012345678901234567890123456789012345678901234567890123456789012
3297: c
3298: c function fluxlim
3299: c
3300: c this function computes the flux limited quick face value
3301: c
3302: c23456789012345678901234567890123456789012345678901234567890123456789012
3305: c#######################################################################
3307: implicit none
3309: double precision fww, fw, fe, fee, vp
3311: double precision fd, fc, fu
3313: double precision f1, f2, f3, f4, fhod, beta, flc
3315: double precision med, quick
3316:
3317: logical limit
3319: c#######################################################################
3321: limit = .true.
3323: if (vp .gt. 0.0d+0) then
3324: fd = fe
3325: fc = fw
3326: fu = fww
3327: else
3328: fd = fw
3329: fc = fe
3330: fu = fee
3331: endif
3333: fhod = quick(fd,fc,fu)
3335: if (limit) then
3337: beta = 0.25d+0
3338: flc = 4.0d+0
3340: f1 = fc
3341: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3342: f3 = fu + ( flc * (fc - fu) )
3343: f4 = med(f1,f2,f3)
3344: fluxlim = vp * med(f1,f4,fhod)
3346: else
3348: fluxlim = vp * fhod
3350: endif
3352: return
3353: end
3354: double precision function fluxlim2(fww,fw,fe,fee,vp)
3355: c23456789012345678901234567890123456789012345678901234567890123456789012
3356: c
3357: c function fluxlim2
3358: c
3359: c this function computes the flux limited quick face value
3360: c
3361: c23456789012345678901234567890123456789012345678901234567890123456789012
3364: c#######################################################################
3366: implicit none
3368: double precision fww, fw, fe, fee, vp
3370: double precision fd, fc, fu
3372: double precision f1, f2, f3, f4, fhod, beta, flc
3374: double precision med, quick
3375:
3376: logical limit, debug
3378: c#######################################################################
3380: debug = .false.
3382: if (debug) then
3383: write(*,*)
3384: write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3385: write(*,*) fee,fe,fw,fww
3386: write(*,*)
3387: endif
3389: limit = .true.
3391: if (vp .gt. 0.0d+0) then
3392: fd = fe
3393: fc = fw
3394: fu = fww
3395: else
3396: fd = fw
3397: fc = fe
3398: fu = fee
3399: endif
3401: fhod = quick(fd,fc,fu)
3403: if (limit) then
3405: beta = 0.25d+0
3406: flc = 4.0d+0
3408: f1 = fc
3409: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3410: f3 = fu + ( flc * (fc - fu) )
3411: f4 = med(f1,f2,f3)
3412: fluxlim2 = med(f1,f4,fhod)
3414: else
3416: fluxlim2 = fhod
3418: endif
3420: return
3421: end
3422: double precision function fprime(x)
3424: implicit none
3426: double precision x, eps
3427: double precision shockp
3429: eps = 1.0d-8
3431: fprime = ( shockp(x+eps) - shockp(x) ) / eps
3433: return
3434: end
3435: double precision function godent(uhl, uhr, hl, hr, eqn)
3436: c23456789012345678901234567890123456789012345678901234567890123456789012
3437: c
3438: c function godent
3439: c
3440: c this function computes the roe/godunov face value plus entropy fix
3441: c
3442: c23456789012345678901234567890123456789012345678901234567890123456789012
3445: c#######################################################################
3447: implicit none
3449: include 'comd.h'
3451: double precision uhl, uhr, hl, hr, ht, uht
3452: double precision lamr1, lamr2, laml1, laml2
3453: double precision deltal1, deltal2
3455: integer eqn
3457: double precision sum
3459: double precision flux
3461: integer i, j
3463:
3465: c#######################################################################
3467: if (debug) then
3468: write(*,*) 'in godent eqn = ',eqn
3469: endif
3471: c do i = 1,neq
3472: c fr(i) = flux(uhr,hr,i)
3473: c fl(i) = flux(uhl,hl,i)
3474: c enddo
3476: deltau(1) = hr - hl
3477: deltau(2) = uhr - uhl
3479: call roestat(uhl, uhr, hl,hr,ht, uht)
3481: call eigen(ht,uht)
3483: do i = 1,neq
3484: sum = 0.0d+0
3485: do j = 1,neq
3486: sum = sum + ( rinv(i,j) * deltau(j) )
3487: enddo
3488: alpha(i) = sum
3489: enddo
3491: deltal1 = 0.0d+0
3492: deltal2 = 0.0d+0
3494: c call eigene(hr,uhr,lamr1, lamr2)
3495: c call eigene(hl,uhl,laml1, laml2)
3496: c
3497: c 1st eigen
3498: c
3499: if ( (laml1 .lt. 0.0d+0) .and.
3500: & (lamr1 .gt. 0.0d+0) ) then
3502: CVAM deltal1 = 4.0d+0 * (lamr1 - laml1)
3503: deltal1 = 4.0d+0 * (lamr1 - laml1) + 1.0d-2
3505: if ( abs(eigval(1)) .lt. (0.5d+0 * deltal1) ) then
3506: eigval(1) = ( ( (eigval(1) ** 2) / deltal1 )
3507: & + ( 0.25d+0 * deltal1 ) )
3508: endif
3510: endif
3511: c
3512: c 2nd eigen
3513: c
3514: if ( (laml2 .lt. 0.0d+0) .and.
3515: & (lamr2 .gt. 0.0d+0) ) then
3517: deltal2 = 4.0d+0 * (lamr2 - laml2)
3519: if ( abs(eigval(2)) .lt. (0.5d+0 * deltal2) ) then
3520: eigval(2) = ( ( (eigval(2) ** 2) / deltal2 )
3521: & + ( 0.25d+0 * deltal2 ) )
3522: endif
3524: endif
3526: if (debug) then
3527: write(*,*)
3528: write(*,*) 'godent debug'
3529: write(*,*) laml1, laml2, lamr1, lamr2
3530: write(*,*) deltal1, deltal2
3531: write(*,*)
3532: endif
3534: do i = 1,neq
3535: sum = 0.0d+0
3536: do j = 1,neq
3537: sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3538: enddo
3539: xnumdif(i) = sum
3540: enddo
3542: do i = 1,neq
3543: froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3544: enddo
3546: godent = froe(eqn)
3548: return
3549: end
3550: double precision function godunov(uhl, uhr, hl, hr, eqn)
3551: c23456789012345678901234567890123456789012345678901234567890123456789012
3552: c
3553: c function godunov
3554: c
3555: c this function computes the roe/godunov face value
3556: c
3557: c23456789012345678901234567890123456789012345678901234567890123456789012
3560: c#######################################################################
3562: implicit none
3564: include 'comd.h'
3566: double precision uhl, uhr, hl, hr, ht, uht
3568: integer eqn
3570: double precision sum
3572: double precision flux
3574: integer i, j
3575:
3577: c#######################################################################
3579: c do i = 1,neq
3580: c fr(i) = flux(uhr,hr,i)
3581: c fl(i) = flux(uhl,hl,i)
3582: c enddo
3584: deltau(1) = hr - hl
3585: deltau(2) = uhr - uhl
3587: call roestat(uhl, uhr, hl,hr,ht, uht)
3589: call eigen(ht,uht)
3591: do i = 1,neq
3592: sum = 0.0d+0
3593: do j = 1,neq
3594: sum = sum + ( rinv(i,j) * deltau(j) )
3595: enddo
3596: alpha(i) = sum
3597: enddo
3599: do i = 1,neq
3600: sum = 0.0d+0
3601: do j = 1,neq
3602: sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3603: enddo
3604: xnumdif(i) = sum
3605: enddo
3607: do i = 1,neq
3608: froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3609: enddo
3611: godunov = froe(eqn)
3614: return
3615: end
3616: double precision function godunov2
3617: & (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3618: c23456789012345678901234567890123456789012345678901234567890123456789012
3619: c
3620: c function godunov2
3621: c
3622: c this function computes the roe/godunov2 face value
3623: c
3624: c23456789012345678901234567890123456789012345678901234567890123456789012
3627: c#######################################################################
3629: implicit none
3631: include 'comd.h'
3632: include 'tube.h'
3634: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3636: integer eqn
3638: double precision rrg, rlg, rurg, rulg, erg, elg
3640: double precision godunov, godent, hlle
3643: c#######################################################################
3645: if (gorder .eq. 1) then
3646: rrg = rr
3647: rlg = rl
3648: rurg = rur
3649: rulg = rul
3650: erg = er
3651: elg = el
3652: else
3653: call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
3654: & rrg, rlg,rurg, rulg, erg, elg)
3655: endif
3657: CVAM if (ientro .eq. 0) then
3658: CVAM godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3659: CVAM elseif(ientro .eq. 1) then
3660: CVAM godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3661: CVAM else
3662: godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3663: CVAM endif
3666: return
3667: end
3668: double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3669: c23456789012345678901234567890123456789012345678901234567890123456789012
3670: c
3671: c function hlle
3672: c
3673: c this function computes the roe/hlle face value
3674: c
3675: c23456789012345678901234567890123456789012345678901234567890123456789012
3678: c#######################################################################
3680: implicit none
3682: include 'comd.h'
3683: include 'tube.h'
3685: double precision rrg,rlg,rurg,rulg,erg,elg
3686: integer eqn
3688: double precision laml1, laml2, laml3
3689: double precision lamr1, lamr2, lamr3
3690: double precision sl, sr
3693: double precision flux
3695: integer i, j, ispeed
3696:
3698: c#######################################################################
3700: ispeed = 1
3702: do i = 1,neq
3703: fr(i) = flux(rrg,rurg,erg,i)
3704: fl(i) = flux(rlg,rulg,elg,i)
3705: enddo
3707: deltau(1) = rrg - rlg
3708: deltau(2) = rurg - rulg
3709: deltau(3) = erg - elg
3711: CVAM call roestat(uhl,uhr,hl,hr,ht,uht)
3713: CVAM call eigene(ht,uht,lamt1, lamt2)
3714: call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3715: call eigene(rlg,rulg,elg,laml1,laml2,laml3)
3717: CVAM if (ispeed .eq. 1) then
3718: CVAM sl = min(laml1,lamt1)
3719: CVAM sr = max(lamt2,lamr2)
3720: CVAM else
3721: sl = min(laml1,lamr1)
3722: sr = max(laml3,lamr3)
3723: CVAM endif
3726: do i = 1,neq
3727: froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )
3728: & / (sr-sl)
3729: enddo
3731: hlle = froe(eqn)
3734: return
3735: end
3736: double precision function med(x1,x2,x3)
3737: c23456789012345678901234567890123456789012345678901234567890123456789012
3738: c
3739: c function med
3740: c
3741: c this function computes the median of three numbers
3742: c
3743: c23456789012345678901234567890123456789012345678901234567890123456789012
3746: c#######################################################################
3748: implicit none
3750: double precision x1, x2, x3
3751: double precision xhi, xlo
3753: c#######################################################################
3755: xhi = max(x1,x2,x3)
3756: xlo = min(x1,x2,x3)
3758: med = x1 + x2 + x3 - xhi - xlo
3760: return
3761: end
3762: double precision function mom
3763: & (rhoee, rhoe, rhop, rhow, rhoww,
3764: & rhouee, rhoue, rhoup, rhouw, rhouww,
3765: & ergee, erge, ergp, ergw, ergww,
3766: & jmom,xold)
3767: c
3768: c This function computes the residual
3769: c for the 1-D momentum equation
3770: c
3771: c
3772: implicit none
3774: include 'comd.h'
3775: include 'tube.h'
3776: c
3777: c input variables
3778: c
3779: double precision rhoee, rhoe, rhop, rhow, rhoww
3780: double precision rhouee, rhoue, rhoup, rhouw, rhouww
3781: double precision ergee, erge, ergp, ergw, ergww
3782: double precision xold(mx*neq)
3783: c
3784: integer jmom
3785: c
3786: c local variables
3787: c
3788: double precision theta1
3789: integer jru
3790: c
3791: c new
3792: c
3793: double precision velfw, velfe
3794: double precision vele,velp,velw
3795: double precision fluxe, fluxw
3796: double precision uurhoe, uurhow
3797: double precision pressee, presse, pressp,pressw, pressww
3798: double precision rupee, rupe, rupp, rupw, rupww
3799: double precision uee, ue, up, uw, uww
3800: double precision source
3801: c
3802: c old
3803: c
3804: double precision rhooee, rhooe, rhoop, rhoow, rhooww
3805: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3806: double precision ergoee, ergoe, ergop, ergow, ergoww
3807: double precision velfow, velfoe
3808: double precision veloe,velop,velow
3809: double precision fluxoe, fluxow
3810: double precision uurhooe, uurhoow
3811: double precision pressoee, pressoe, pressop, pressow, pressoww
3812: double precision rupoee, rupoe, rupop, rupow, rupoww
3813: double precision uoee, uoe, uop, uow, uoww
3814: double precision sourceo
3816: double precision eps
3817: c
3818: c functions
3819: c
3820: double precision godunov2, eos
3821: double precision upwind, fluxlim
3822: c
3823: c
3824: c ******************************************************************
3825: c
3826: c
3827: eps = 1.0d-32
3828: c
3829: jru = (neq*jmom) - 1
3831: c########################
3832: c
3833: c NEW
3834: c
3835: c########################
3837: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
3838: & rhouee, rhoue, rhoup, rhouw, rhouww,
3839: & ergee, erge, ergp, ergw, ergww,
3840: & vele, velp, velw, jmom)
3842: presse = eos(rhoe, rhoue, erge )
3843: pressw = eos(rhow, rhouw, ergw )
3845: velfe = 0.5d+0 * (vele + velp)
3846: velfw = 0.5d+0 * (velw + velp)
3848: if (ihod .eq. 1) then
3850: uurhoe = upwind(rhoup,rhoue,velfe)
3851: uurhow = upwind(rhouw,rhoup,velfw)
3853: elseif (ihod .eq. 2) then
3855: uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3856: uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)
3858: endif
3860: if (ihod .eq. 3) then
3861: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
3862: & rhouw,rhoup,rhoue,rhouee,
3863: & ergw, ergp, erge, ergee,2)
3864: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
3865: & rhouww,rhouw,rhoup,rhoue,
3866: & ergww, ergw, ergp, erge,2)
3867: else
3868: fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3869: fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3870: endif
3873: source = 0.0d+0
3875: c########################
3876: c
3877: c OLD
3878: c
3879: c########################
3881: call Setpbc(jmom,xold,
3882: & rhooee, rhooe, rhoop, rhoow, rhooww,
3883: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
3884: & ergoee, ergoe, ergop, ergow, ergoww,
3885: & veloe, velop, velow)
3887: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
3888: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
3889: & ergoee, ergoe, ergop, ergow, ergoww,
3890: & veloe, velop, velow, jmom)
3892: pressoe = eos(rhooe, rhouoe, ergoe)
3893: pressow = eos(rhoow, rhouow, ergow)
3895: velfoe = 0.5d+0 * (veloe + velop)
3896: velfow = 0.5d+0 * (velow + velop)
3898: if (ihod .eq. 1) then
3900: uurhooe = upwind(rhouop,rhouoe,velfoe)
3901: uurhoow = upwind(rhouow,rhouop,velfow)
3903: elseif (ihod .eq. 2) then
3905: uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3906: uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)
3908: endif
3910: if (ihod .eq. 3) then
3911: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3912: & rhouow,rhouop,rhouoe,rhouoee,
3913: & ergow, ergop, ergoe, ergoee,2)
3914: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3915: & rhouoww,rhouow,rhouop,rhouoe,
3916: & ergoww, ergow, ergop, ergoe,2)
3917: else
3918: fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3919: fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3920: endif
3922: sourceo = 0.0d+0
3925: c########################
3926: c
3927: c FUNCTION
3928: c
3929: c########################
3931: theta1 = 1.0d+0 - theta
3932: mom = (rhoup - xold(jru))
3933: & + ( theta * ( (fluxe - fluxw ) - source ) )
3934: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3935: CVAM
3936: if (probnum .eq. 3) then
3937: mom = 0.0d+0
3938: endif
3939: CVAM
3940: if (debug) then
3941: write(*,*)
3942: write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3943: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3944: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3945: write(*,*) 'presse = ',presse,'pressw = ',pressw
3946: write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3947: write(*,*)
3948: endif
3950: return
3951: end
3952: double precision function quick(fd, fc, fu)
3953: c23456789012345678901234567890123456789012345678901234567890123456789012
3954: c
3955: c function quick
3956: c
3957: c this function computes the quick face value
3958: c
3959: c23456789012345678901234567890123456789012345678901234567890123456789012
3962: c#######################################################################
3964: implicit none
3966: double precision fd, fc, fu
3968: c#######################################################################
3970: quick = ( (3.0d+0 * fd) + (6.0d+0 * fc) - fu ) / 8.0d+0
3972: return
3973: end
3974: double precision function rexact(x,t)
3976: implicit none
3978: double precision x,t
3979: double precision xot, head, tail, contact, ufan
3980: double precision xpow, grat, urat
3981: double precision uexact
3984: logical debug
3986: include 'tube.h'
3988: debug = .false.
3991: if (t .le. 0.0d+0) then
3992: if (x .gt. 0.0d+0) then
3993: rexact = r1
3994: else
3995: rexact = r4
3996: endif
3997: else
3999: xot = x/t
4000: head = -a4
4001: tail = v3 - a3
4002: contact = v2
4004: if (xot .lt. head) then
4005: rexact = r4
4006: elseif (xot .gt. sspd) then
4007: rexact = r1
4008: elseif (xot .gt. contact) then
4009: rexact = r2
4010: elseif (xot .gt. tail) then
4011: rexact = r3
4012: else
4013: ufan = uexact(x,t)
4014: grat = (gamma - 1.0d+0) / 2.0d+0
4015: xpow = 1.0d+0 / grat
4016: urat = ufan / a4
4017: rexact = r4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
4018: endif
4020: endif
4023: if (debug) then
4024: write(*,*)
4025: write(*,*) 'rexact(',x,',',t,') = ',rexact
4026: write(*,*)
4027: endif
4029: return
4030: end
4031: subroutine roestat(uhl, uhr, hl,hr,ht,uht)
4032: c23456789012345678901234567890123456789012345678901234567890123456789012
4033: c
4034: c subroutine roestat
4035: c
4036: c This subroutine computes the roe state at a face
4037: c
4038: c23456789012345678901234567890123456789012345678901234567890123456789012
4041: c#######################################################################
4043: implicit none
4045: include 'comd.h'
4047: double precision uhl, uhr, hl, hr, ht, uht
4049: double precision ul, ur, shl, shr, xnum, xdenom
4050:
4053: c#######################################################################
4055: ul = uhl / hl
4056: ur = uhr / hr
4058: shl = sqrt(hl)
4059: shr = sqrt(hr)
4061: xnum = (shl * ul) + (shr * ur)
4062: xdenom = shl + shr
4064: ht = 0.5d+0 * (hl + hr)
4065: uht = ht * ( xnum / xdenom )
4067: return
4068: end
4069: subroutine rval2
4071: implicit none
4073: double precision prat, grat, xnum, xdenom
4076: logical debug
4078: include 'tube.h'
4080: debug = .false.
4082: prat = p2/p1
4083: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
4085: xnum = 1.0d+0 + (grat * prat)
4086: xdenom = grat + prat
4087:
4088: r2 = r1 * (xnum/xdenom)
4089:
4092: if (debug) then
4093: write(*,*)
4094: write(*,*) 'r1 = ',r1
4095: write(*,*) 'r2 = ',r2
4096: endif
4098: return
4099: end
4100: subroutine secondq
4101: & (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
4102: & rrg, rlg,rurg, rulg, erg, elg)
4103: c23456789012345678901234567890123456789012345678901234567890123456789012
4104: c
4105: c subroutine secondq
4106: c
4107: c this subroutine computes the second order (based on quick) left
4108: c and right states for the godunov solver.
4109: c
4110: c23456789012345678901234567890123456789012345678901234567890123456789012
4113: c#######################################################################
4115: implicit none
4117: include 'comd.h'
4119: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
4120: double precision rrg, rlg,rurg, rulg, erg, elg
4124: double precision veld, ull,ul,ur,urr, ulg, urg
4126: double precision fluxlim2
4129: c#######################################################################
4131: c
4132: c compute the velocities
4133: c
4134: ull = rull/rll
4135: ul = rul /rl
4136: ur = rur /rr
4137: urr = rurr/rrr
4139: c
4140: c compute the left state first
4141: c
4142: veld = 1.0d+0
4144: rlg = fluxlim2(rll,rl,rr,rrr,veld)
4145: ulg = fluxlim2(ull,ul,ur,urr,veld)
4146: rulg = rlg * ulg
4147: elg = fluxlim2(ell,el,er,err,veld)
4148: c
4149: c now compute the right state
4150: c
4151: veld = -1.0d+0
4153: rrg = fluxlim2(rll,rl,rr,rrr,veld)
4154: urg = fluxlim2(ull,ul,ur,urr,veld)
4155: rurg = rrg * urg
4156: erg = fluxlim2(ell,el,er,err,veld)
4160: return
4161: end
4162: double precision function shockp(x)
4164: implicit none
4166: double precision x
4167: double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
4168: logical debug
4170: include 'tube.h'
4172: debug = .false.
4175: if (debug) then
4176: write(*,*)
4177: write(*,*) 'gamma = ',gamma
4178: write(*,*) 'a1 = ',a1
4179: write(*,*) 'a4 = ',a4
4180: write(*,*) 'p1 = ',p1
4181: write(*,*) 'p2 = ',x
4182: write(*,*)
4183: endif
4185: xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
4186: xdenom = sqrt ( 2.0d+0 * gamma * ( (2.0d+0*gamma)
4187: & + (gamma + 1.0d+0) * ((x/p1) - 1) ) )
4188: xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)
4190: shockp = (x/p1)*((1.0d+0-(xnum/xdenom))**xpow) - (p4/p1)
4193: if (debug) then
4194: write(*,*)
4195: write(*,*) 'xnum = ',xnum
4196: write(*,*) 'gamma = ',gamma
4197: write(*,*) 'a1 = ',a1
4198: write(*,*) 'a4 = ',a4
4199: write(*,*) 'p1 = ',p1
4200: write(*,*) 'xdenom = ',xdenom
4201: write(*,*) 'xpow = ',xpow
4202: write(*,*) 'shockp = ',shockp
4203: write(*,*) 'p2 = ',x
4204: write(*,*)
4205: endif
4207: return
4208: end
4209: double precision function uexact(x,t)
4211: implicit none
4213: double precision x,t
4214: double precision xot, head, tail
4217: logical debug
4219: include 'tube.h'
4221: debug = .false.
4223: if (debug) then
4224: write(*,*)
4225: write(*,*) 't = ',t
4226: write(*,*) 'x = ',x
4227: write(*,*) 'a4 = ',a4
4228: write(*,*) 'v3 = ',v3
4229: write(*,*) 'a3 = ',a3
4230: write(*,*)
4231: endif
4233: if (t .le. 0.0d+0) then
4234: uexact = 0.0d+0
4235: else
4237: xot = x/t
4238: head = -a4
4239: tail = v3 - a3
4241: if (xot .lt. head) then
4242: uexact = 0.0d+0
4243: elseif (xot .gt. sspd) then
4244: uexact = 0.0d+0
4245: elseif (xot .gt. tail) then
4246: uexact = v2
4247: else
4248: uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
4249: endif
4251: endif
4254: if (debug) then
4255: write(*,*)
4256: CVAM write(*,*) 'x = ',x,' t = ',t
4257: write(*,*) 'uexact = ',uexact
4258: write(*,*)
4259: endif
4261: return
4262: end
4263: double precision function upwind(fw, fe, vp)
4264: c23456789012345678901234567890123456789012345678901234567890123456789012
4265: c
4266: c function upwind
4267: c
4268: c this function computes the upwind face value
4269: c
4270: c23456789012345678901234567890123456789012345678901234567890123456789012
4273: c#######################################################################
4275: implicit none
4277: double precision fw, fe, vp
4279: c#######################################################################
4281: if (vp .gt. 0.0) then
4282: upwind = vp * fw
4283: else
4284: upwind = vp * fe
4285: endif
4287: return
4288: end
4289: subroutine uval2
4291: implicit none
4293: double precision prat, grat1, grat2, arat, xnum
4296: logical debug
4298: include 'tube.h'
4300: debug = .false.
4302: prat = p2/p1
4303: grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
4304: grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
4305: arat = a1/gamma
4307: xnum = sqrt ( grat2 / (prat + grat1) )
4309: v2 = arat * (prat - 1.0d+0) * xnum
4311: if (debug) then
4312: write(*,*)
4313: write(*,*) 'v2 = ',v2
4314: endif
4316: return
4317: end
4318: subroutine val3
4320: implicit none
4322: double precision prat, rpow, epow, p3t
4325: logical debug
4327: include 'tube.h'
4329: debug = .false.
4332: p3 = p2
4334: prat = p3/p4
4336: rpow = 1.0d+0 / gamma
4338: r3 = r4 * ( prat ** rpow )
4340: epow = (gamma - 1.0d+0) / gamma
4342: e3 = e4 * ( (p3/p4) ** epow )
4344: p3t = (gamma - 1.0d+0) * r3 * e3
4346: a3 = sqrt(gamma*p3/r3)
4348: if (debug) then
4349: write(*,*)
4350: write(*,*) 'a3 = ',a3
4351: write(*,*) 'r3 = ',r3
4352: write(*,*) 'e3 = ',e3
4353: write(*,*) 'p3 = ',p3
4354: write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
4355: write(*,*)
4356: endif
4358: return
4359: end
4360: subroutine wval
4362: implicit none
4364: double precision prat, grat, xnum
4367: logical debug
4369: include 'tube.h'
4371: debug = .false.
4373: prat = p2/p1
4374: grat = (gamma + 1.0d+0) / (2.0d+0 * gamma)
4376: xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
4377:
4378: sspd = a1 * sqrt(xnum)
4379:
4382: if (debug) then
4383: write(*,*)
4384: write(*,*) 'sspd = ',sspd
4385: endif
4387: return
4388: end