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: ! petsc.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/petsc.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/petsc.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/petsc.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/petsc.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/petsc.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*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/petsc.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/petsc.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/petsc.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(ctx,x,y,ierr)
1888: ! -------------------------------------------------------------------
1889: !
1890: ! PCRadApply - This routine demonstrates the use of a
1891: ! user-provided preconditioner.
1892: !
1893: ! Input Parameters:
1894: ! dummy - optional user-defined context, not used here
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/petsc.h
1916: #include finclude/petscvec.h
1918: #include "comd.h"
1919: c
1920: c Input
1921: c
1922: PetscFortranAddr ctx(*)
1923: Vec x, y
1924: integer ierr
1925: c
1926: c Local
1927: c
1928: IS defrom, deto
1929: Vec de, rese
1930: VecScatter scatter
1931: PetscScalar lde_v(0:1),lrese_v(0:1)
1932: PetscOffset lde_i, lrese_i
1933: c
1934: c Identity preconditioner
1935: c
1936: call VecCopy(x,y,ierr)
1937: c
1938: c if kappa0 not equal to zero then precondition the radiation diffusion
1939: c
1940: if (kappa0 .ne. 0.0d+0) then
1941:
1943: c
1944: c Create needed vectors
1945: c
1946: CALL VecCreate(PETSC_COMM_WORLD,de,ierr)
1947: CALL VecSetSizes(de,PETSC_DECIDE,mx,ierr)
1948: CALL VecSetType(de,VECMPI,ierr)
1949: call VecSetFromOptions(de,ierr)
1950: call VecDuplicate(de,rese,ierr)
1951: c
1952: c create index sets for scatters
1953: c
1954: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1955: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1956: c
1957: c load rese from x
1958: c
1959: call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1960: call VecScatterBegin(scatter,x,rese,INSERT_VALUES,
1961: & SCATTER_FORWARD,ierr)
1962: call VecScatterEnd(scatter,x,rese,INSERT_VALUES,
1963: & SCATTER_FORWARD,ierr)
1964: call VecScatterDestroy(scatter,ierr)
1965: c
1966: c apply preconditioner
1967: c
1968: call PCApply(ctx(1),rese,de,ierr)
1970: if (debug) then
1971: write(*,*) 'PCRadApply dh is'
1972: call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1973: endif
1974: c
1975: c load de into y
1976: c
1977: call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1978: call VecScatterBegin(scatter,de,y,INSERT_VALUES,
1979: & SCATTER_FORWARD,ierr)
1980: call VecScatterEnd(scatter,de,y,INSERT_VALUES,
1981: & SCATTER_FORWARD,ierr)
1982: call VecScatterDestroy(scatter,ierr)
1984: if (debug) then
1985: write(*,*) 'PCRadApply y is'
1986: call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1987: endif
1989: call VecDestroy(de,ierr)
1990: call VecDestroy(rese,ierr)
1992: call ISDestroy(defrom,ierr)
1993: call ISDestroy(deto,ierr)
1995: endif
1998: return
1999: end
2000: subroutine PCRadSetUp(ctx,ierr)
2001: !
2002: ! PCRadSetUp - This routine sets up a user-defined
2003: ! preconditioner context.
2004: !
2005: ! Input Parameters:
2006: ! In this example the shell preconditioner application context
2007: ! is a Fortran integer array:
2008: ! ctx(1) = shell preconditioner pressure matrix contex
2009: ! ctx(2) = semi implicit pressure matrix
2010: ! ctx(4) = xold - old time values need for time advancement
2011: ! ctx(5) = mx - number of control volumes
2012: ! ctx(6) = N - total number of unknowns
2013: !
2014: ! Output Parameter:
2015: ! ierr - error code (nonzero if error has been detected)
2016: !
2017: ! Notes:
2018: ! In this example, we define the shell preconditioner to be Jacobi
2019: ! method. Thus, here we create a work vector for storing the reciprocal
2020: ! of the diagonal of the preconditioner matrix; this vector is then
2021: ! used within the routine PCRadApply().
2022: !
2024: implicit none
2026: #include finclude/petsc.h
2027: #include finclude/petscvec.h
2028: #include finclude/petscmat.h
2030: #include "comd.h"
2031: c
2032: c Input
2033: c
2034: PetscFortranAddr ctx(*)
2035: integer ierr
2036: c
2037: c Local
2038: c
2039: Vec eold
2040:
2041: PetscScalar le_v(0:1)
2042: PetscOffset le_i
2043:
2044: c
2045: c create vector
2046: c
2047: CALL VecCreate(PETSC_COMM_WORLD,eold,ierr)
2048: CALL VecSetSizes(eold,PETSC_DECIDE,mx,ierr)
2049: CALL VecSetType(eold,VECMPI,ierr)
2050: call VecSetFromOptions(eold,ierr)
2051: c
2052: c set up the matrix based on xold
2053: c
2054: call Setmat(ctx,ierr)
2055: c
2056: c set up the preconditioner
2057: c
2058: call PCDestroy(ctx(1),ierr)
2059: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
2060: CVAM call PCSetType(ctx(1),PCJACOBI,ierr)
2061: call PCSetType(ctx(1),PCLU,ierr)
2062: ! call PCSetVector(ctx(1),eold,ierr)
2063: call PCSetOperators(ctx(1),ctx(2),ctx(2),
2064: & DIFFERENT_NONZERO_PATTERN,ierr)
2065: call PCSetUp(ctx(1),ierr)
2067: call VecDestroy(eold,ierr)
2070: return
2071: end
2072: subroutine Setmat(ctx,ierr)
2074: implicit none
2076: #include finclude/petsc.h
2077: #include finclude/petscvec.h
2078: #include finclude/petscmat.h
2079: ! Common blocks:
2080: #include "comd.h"
2081: #include "tube.h"
2083: ! Input/output variables:
2084: PetscFortranAddr ctx(*)
2085: integer ierr
2087: ! Local variables:
2088: PetscScalar lx_v(0:1)
2089: PetscOffset lx_i
2091: double precision xmult, himh, hiph, diag, upper, lower
2092: double precision hi, hip1, him1
2093: double precision
2094: & rhoee, rhoe, rhop, rhow, rhoww,
2095: & rhouee, rhoue, rhoup, rhouw, rhouww,
2096: & ergee, erge, ergp, ergw, ergww,
2097: & ue, up, uw
2098: double precision see, sep, sew, seef, sewf, tef, twf,
2099: & ref, rwf, kef, kwf, xmulte, xmultw
2100: c
2101: integer im, nx, ny
2102: c
2103: c get pointers to xold
2104: c
2105: call VecGetArray(ctx(4),lx_v,lx_i,ierr)
2106:
2108: c
2109: c############################
2110: c
2111: c loop over all cells begin
2112: c
2113: c############################
2114: c
2115: do im = 1,mx
2116: c
2117: c set scalars
2118: c
2119: call Setpbc(im,lx_v(lx_i),
2120: & rhoee, rhoe, rhop, rhow, rhoww,
2121: & rhouee, rhoue, rhoup, rhouw, rhouww,
2122: & ergee, erge, ergp, ergw, ergww,
2123: & ue, up, uw)
2124: c
2125: c set diffusion coefficients
2126: c
2127: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2128: sep = (ergp/rhop) - (0.5d+0 * up * up)
2129: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2131: seef = 0.5d+0 * (see + sep)
2132: sewf = 0.5d+0 * (sew + sep)
2134: tef = seef / csubv
2135: twf = sewf / csubv
2137: ref = 0.5d+0 * (rhoe + rhop)
2138: rwf = 0.5d+0 * (rhow + rhop)
2140: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2141: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2143: if (wilson) then
2144: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2145: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2146: endif
2147: c
2148: c set coefficients
2149: c
2150: xmult = dt / (dx * dx * csubv)
2152: xmulte = xmult * kef
2153: xmultw = xmult * kwf
2155: upper = -(xmulte / rhoe)
2156: lower = -(xmultw / rhow)
2158: diag = 1.0d+0 + ( (xmulte + xmultw) / rhop )
2160: c
2161: c load coefficients into the matrix
2162: c
2163: call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)
2165: if (im .eq. 1) then
2166: call MatSetValues(ctx(2),1,im-1,1,im ,upper,
2167: 1 INSERT_VALUES,ierr)
2168: elseif (im .eq. mx) then
2169: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2170: 1 INSERT_VALUES,ierr)
2171: else
2172: call MatSetValues(ctx(2),1,im-1,1,im ,upper,
2173: 1 INSERT_VALUES,ierr)
2174: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,
2175: 1 INSERT_VALUES,ierr)
2176: endif
2179: enddo
2180: c
2181: c############################
2182: c
2183: c loop over all cells end
2184: c
2185: c############################
2186: c
2187:
2188: c
2189: c final load of matrix
2190: c
2191: call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2192: call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2194: if (debug) then
2195: call MatGetSize(ctx(2),nx,ny,ierr)
2196: write(*,*) 'in setup nx = ',nx,' ny = ',ny
2197: call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2198: endif
2200: call VecRestoreArray (ctx(4),lx_v,lx_i,ierr)
2204: return
2205: end
2206: subroutine Setpbc(i,x,
2207: & rhoee, rhoe, rhop, rhow, rhoww,
2208: & rhouee, rhoue, rhoup, rhouw, rhouww,
2209: & ergee, erge, ergp, ergw, ergww,
2210: & vele, velp, velw)
2212: implicit none
2214: ! Common blocks:
2215: #include "comd.h"
2217: ! Input/output variables:
2218: PetscScalar x(mx*neq)
2219: integer i
2220: double precision rhoee, rhoe, rhop, rhow, rhoww
2221: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2222: double precision ergee, erge, ergp, ergw, ergww
2223: double precision vele, velp, velw
2225: ! Local variables:
2226: integer jr, jru, je
2228: c
2229: c set pointers
2230: c
2231: jr = (neq*i) - 2
2232: jru = (neq*i) - 1
2233: je = (neq*i)
2235: if (debug) then
2236: write(*,*)
2237: write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2238: write(*,*)
2239: endif
2240:
2241: if (i .eq. 1) then
2243: rhoee = x(jr+(2*neq))
2244: rhoe = x(jr+neq)
2245: rhop = x(jr)
2246: rhow = x(jr)
2247: rhoww = x(jr)
2249: rhouee = x(jru+(2*neq))
2250: rhoue = x(jru+neq)
2251: rhoup = x(jru)
2252: rhouw = x(jru)
2253: rhouww = x(jru)
2255: ergee = x(je+(2*neq))
2256: erge = x(je+neq)
2257: ergp = x(je)
2258: ergw = x(je)
2259: ergww = x(je)
2261: velw = 0.0d+0
2262: velp = rhoup/rhop
2263: vele = rhoue/rhoe
2265: elseif (i .eq. 2) then
2267: rhoee = x(jr+(2*neq))
2268: rhoe = x(jr+neq)
2269: rhop = x(jr)
2270: rhow = x(jr-neq)
2271: rhoww = x(jr-neq)
2273: rhouee = x(jru+(2*neq))
2274: rhoue = x(jru+neq)
2275: rhoup = x(jru)
2276: rhouw = x(jru-neq)
2277: rhouww = x(jru-neq)
2279: ergee = x(je+(2*neq))
2280: erge = x(je+neq)
2281: ergp = x(je)
2282: ergw = x(je-neq)
2283: ergww = x(je-neq)
2285: velw = rhouw/rhow
2286: velp = rhoup/rhop
2287: vele = rhoue/rhoe
2289: elseif (i .eq. mx-1) then
2291: rhoee = x(jr+neq)
2292: rhoe = x(jr+neq)
2293: rhop = x(jr)
2294: rhow = x(jr-neq)
2295: rhoww = x(jr-(2*neq))
2297: rhouee = x(jru+neq)
2298: rhoue = x(jru+neq)
2299: rhoup = x(jru)
2300: rhouw = x(jru-neq)
2301: rhouww = x(jru-(2*neq))
2303: ergee = x(je+neq)
2304: erge = x(je+neq)
2305: ergp = x(je)
2306: ergw = x(je-neq)
2307: ergww = x(je-(2*neq))
2309: velw = rhouw/rhow
2310: velp = rhoup/rhop
2311: vele = rhoue/rhoe
2313: elseif (i .eq. mx) then
2315: rhoee = x(jr)
2316: rhoe = x(jr)
2317: rhop = x(jr)
2318: rhow = x(jr-neq)
2319: rhoww = x(jr-(2*neq))
2321: rhouee = x(jru)
2322: rhoue = x(jru)
2323: rhoup = x(jru)
2324: rhouw = x(jru-neq)
2325: rhouww = x(jru-(2*neq))
2327: ergee = x(je)
2328: erge = x(je)
2329: ergp = x(je)
2330: ergw = x(je-neq)
2331: ergww = x(je-(2*neq))
2333: velw = rhouw/rhow
2334: velp = rhoup/rhop
2335: vele = 0.0d+0
2337: else
2339: rhoee = x(jr+(2*neq))
2340: rhoe = x(jr+neq)
2341: rhop = x(jr)
2342: rhow = x(jr-neq)
2343: rhoww = x(jr-(2*neq))
2345: rhouee = x(jru+(2*neq))
2346: rhoue = x(jru+neq)
2347: rhoup = x(jru)
2348: rhouw = x(jru-neq)
2349: rhouww = x(jru-(2*neq))
2351: ergee = x(je+(2*neq))
2352: erge = x(je+neq)
2353: ergp = x(je)
2354: ergw = x(je-neq)
2355: ergww = x(je-(2*neq))
2357: velw = rhouw/rhow
2358: velp = rhoup/rhop
2359: vele = rhoue/rhoe
2361: endif
2363: if (debug) then
2364: write(*,*)
2365: write(*,*) 'in Setpbc ',i,jr,jru,je
2366: write(*,*) 'mx = ',mx
2367: write(*,*)
2368: endif
2371: return
2372: end
2373: subroutine Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2374: & rhouee, rhoue, rhoup, rhouw, rhouww,
2375: & ergee, erge, ergp, ergw, ergww,
2376: & ue, up, uw, jbc)
2378: implicit none
2380: ! Common blocks:
2381: #include "comd.h"
2383: ! Input/output variables:
2384: integer jbc
2385: double precision rhoee, rhoe, rhop, rhow, rhoww
2386: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2387: double precision ergee, erge, ergp, ergw, ergww
2388: double precision ue, up, uw
2390: ! Local variables:
2392: if (jbc .eq. 1) then
2393: rhoww = rhop
2394: rhow = rhop
2395: rhouww = rhoup
2396: rhouw = rhoup
2397: ergww = ergp
2398: ergw = ergp
2399: uw = 0.0d+0
2400: elseif (jbc .eq. 2) then
2401: rhoww = rhow
2402: rhouww = rhouw
2403: ergww = ergw
2404: uw = rhouw / rhow
2405: else
2406: uw = rhouw / rhow
2407: endif
2409: if (jbc .eq. mx) then
2410: rhoee = rhop
2411: rhoe = rhop
2412: rhouee = rhoup
2413: rhoue = rhoup
2414: ergee = ergp
2415: erge = ergp
2416: ue = 0.0d+0
2417: elseif (jbc .eq. mx-1) then
2418: rhoee = rhoe
2419: rhouee = rhoue
2420: ergee = erge
2421: ue = rhoue / rhoe
2422: else
2423: ue = rhoue / rhoe
2424: endif
2426: up = rhoup / rhop
2428: if (debug) then
2429: write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2430: endif
2433: return
2434: end
2435: double precision function cont
2436: & (rhoee, rhoe, rhop, rhow, rhoww,
2437: & rhouee, rhoue, rhoup, rhouw, rhouww,
2438: & ergee, erge, ergp, ergw, ergww,
2439: & jcont,xold)
2440: c
2441: c This function computes the residual
2442: c for the 1-D continuity equation
2443: c
2444: c
2445: implicit none
2447: include 'comd.h'
2448: include 'tube.h'
2449: c
2450: c input variables
2451: c
2452: double precision rhoee, rhoe, rhop, rhow, rhoww
2453: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2454: double precision ergee, erge, ergp, ergw, ergww
2455: double precision xold(mx*neq)
2456: c
2457: integer jcont
2458: c
2459: c local variables
2460: c
2461: double precision theta1
2462: integer jr
2463: c
2464: c new
2465: c
2466: double precision velfw, velfe
2467: double precision vele,velp,velw
2468: double precision fluxe, fluxw
2469: double precision urhoe, urhow
2470: double precision source
2471: c
2472: c old
2473: c
2474: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2475: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2476: double precision ergoee, ergoe, ergop, ergow, ergoww
2477: double precision teoee, teoe, teop, teow, teoww,
2478: & uoe, uoee, uop, uow, uoww
2479: double precision velfow, velfoe
2480: double precision veloe,velop,velow
2481: double precision fluxoe, fluxow
2482: double precision urhooe, urhoow
2483: double precision sourceo
2484: c
2485: c functions
2486: c
2487: double precision godunov2
2488: double precision upwind, fluxlim
2489: c
2490: c
2491: c ******************************************************************
2492: c
2493: c
2494: c
2495: if (debug) then
2496: write(*,*)
2497: write(*,*) 'in cont',jcont,' ihod = ',ihod
2498: write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2499: write(*,*) 'rhop = ',rhop
2500: write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2501: write(*,*)
2502: endif
2504: jr = (neq*jcont) - 2
2506: c########################
2507: c
2508: c NEW
2509: c
2510: c########################
2512: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2513: & rhouee, rhoue, rhoup, rhouw, rhouww,
2514: & ergee, erge, ergp, ergw, ergww,
2515: & vele, velp, velw, jcont)
2517: velfe = 0.5d+0 * (vele + velp)
2518: velfw = 0.5d+0 * (velw + velp)
2520: if (ihod .eq. 1) then
2522: urhoe = upwind(rhop,rhoe,velfe)
2523: urhow = upwind(rhow,rhop,velfw)
2525: elseif (ihod .eq. 2) then
2527: urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2528: urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)
2530: endif
2532: if (ihod .eq. 3) then
2533: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2534: & rhouw,rhoup,rhoue,rhouee,
2535: & ergw, ergp, erge, ergee,1)
2536: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2537: & rhouww,rhouw,rhoup,rhoue,
2538: & ergww, ergw, ergp, erge,1)
2539: else
2540: fluxe = (dt/dx) * urhoe
2541: fluxw = (dt/dx) * urhow
2542: endif
2545: source = 0.0d+0
2547: c########################
2548: c
2549: c OLD
2550: c
2551: c########################
2553: call Setpbc(jcont,xold,
2554: & rhooee, rhooe, rhoop, rhoow, rhooww,
2555: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2556: & ergoee, ergoe, ergop, ergow, ergoww,
2557: & veloe, velop, velow)
2559: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
2560: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2561: & ergoee, ergoe, ergop, ergow, ergoww,
2562: & veloe, velop, velow, jcont)
2564: velfoe = 0.5d+0 * (veloe + velop)
2565: velfow = 0.5d+0 * (velow + velop)
2568: if (ihod .eq. 1) then
2570: urhooe = upwind(rhoop,rhooe,velfoe)
2571: urhoow = upwind(rhoow,rhoop,velfow)
2573: elseif (ihod .eq. 2) then
2575: urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2576: urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)
2578: endif
2580: if (ihod .eq. 3) then
2581: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
2582: & rhouow,rhouop,rhouoe,rhouoee,
2583: & ergow, ergop, ergoe, ergoee,1)
2584: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
2585: & rhouoww,rhouow,rhouop,rhouoe,
2586: & ergoww, ergow, ergop, ergoe,1)
2587: else
2588: fluxoe = (dt/dx) * urhooe
2589: fluxow = (dt/dx) * urhoow
2590: endif
2592: sourceo = 0.0d+0
2595: c########################
2596: c
2597: c FUNCTION
2598: c
2599: c########################
2601: theta1 = 1.0d+0 - theta
2602: cont = (rhop - xold(jr))
2603: & + ( theta * ( (fluxe - fluxw ) - source ) )
2604: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2605: CVAM
2606: if (probnum .eq. 3) then
2607: cont = 0.0d+0
2608: endif
2609: CVAM
2612: if (debug) then
2613: write(*,*)
2614: write(*,*) 'cont(',jcont,') = ',cont
2615: write(*,*) 'theta = ',theta,'rhop = ',rhop
2616: write(*,*) 'source = ',source,' sourceo = ',sourceo
2617: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2618: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2619: write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2620: write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2621: write(*,*)
2622: endif
2624: return
2625: end
2626: double precision function eexact(x,t)
2628: implicit none
2630: double precision x,t
2631: double precision xot, head, tail, contact, ufan
2632: double precision xpow, grat, urat
2633: double precision uexact
2636: logical debug
2638: include 'tube.h'
2640: debug = .false.
2643: if (t .le. 0.0d+0) then
2644: if (x .gt. 0.0d+0) then
2645: eexact = e1
2646: else
2647: eexact = e4
2648: endif
2649: else
2651: xot = x/t
2652: head = -a4
2653: tail = v3 - a3
2654: contact = v2
2656: if (xot .lt. head) then
2657: eexact = e4
2658: elseif (xot .gt. sspd) then
2659: eexact = e1
2660: elseif (xot .gt. contact) then
2661: eexact = e2
2662: elseif (xot .gt. tail) then
2663: eexact = e3
2664: else
2665: ufan = uexact(x,t)
2666: grat = (gamma - 1.0d+0) / 2.0d+0
2667: xpow = 2.0d+0
2668: urat = ufan / a4
2669: eexact = e4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
2670: endif
2672: endif
2675: if (debug) then
2676: write(*,*)
2677: write(*,*) 'eexact(',x,',',t,') = ',eexact
2678: write(*,*)
2679: endif
2681: return
2682: end
2683: subroutine eigen(ht,uht)
2684: c23456789012345678901234567890123456789012345678901234567890123456789012
2685: c
2686: c subroutine eigen
2687: c
2688: c This subroutine computes the eigen values and eigen vectors
2689: c
2690: c23456789012345678901234567890123456789012345678901234567890123456789012
2693: c#######################################################################
2695: implicit none
2697: include 'comd.h'
2699: double precision ht, uht
2701: double precision ut, at, lam1, lam2
2704: c#######################################################################
2706: ut = uht / ht
2707: at = sqrt( ht)
2709: lam1 = ut - at
2710: lam2 = ut + at
2712: eigval(1) = lam1
2713: eigval(2) = lam2
2715: eigvec(1,1) = 1.0d+0
2716: eigvec(2,1) = lam1
2717: eigvec(1,2) = 1.0d+0
2718: eigvec(2,2) = lam2
2720: rinv(1,1) = lam2 / (2.0d+0 * at)
2721: rinv(2,1) = -lam1 / (2.0d+0 * at)
2722: rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2723: rinv(2,2) = 1.0d+0 / (2.0d+0 * at)
2726: return
2727: end
2728: subroutine eigene(r,ru,e,l1,l2,l3)
2729: c23456789012345678901234567890123456789012345678901234567890123456789012
2730: c
2731: c subroutine eigene
2732: c
2733: c This subroutine computes the eigen values for the entropy fix
2734: c
2735: c23456789012345678901234567890123456789012345678901234567890123456789012
2738: c#######################################################################
2740: implicit none
2742: include 'tube.h'
2744: double precision r,ru,e,l1,l2,l3
2746: double precision p,u,a
2748: double precision eos
2749: integer ierr
2751: logical debug
2754: c#######################################################################
2756: debug = .false.
2758: if (debug) then
2759: write(*,*)
2760: write(*,*) 'gamma = ',gamma
2761: write(*,*) 'r,ru,e = ',r,ru,e
2762: write(*,*)
2763: endif
2765: p = eos(r,ru,e)
2766: u = ru/r
2767: if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2768: write(*,*)
2769: write(*,*) 'gamma = ',gamma
2770: write(*,*) 'r = ',r
2771: write(*,*) 'p = ',p
2772: write(*,*)
2773: call PetscFinalize(ierr)
2774: stop
2775: endif
2776: a = sqrt((gamma * p)/r)
2778: if (debug) then
2779: write(*,*)
2780: write(*,*) 'p,u,a = ',p,u,a
2781: write(*,*)
2782: endif
2784: l1 = u - a
2785: l2 = u
2786: l3 = u + a
2788: return
2789: end
2790: double precision function energy
2791: & (rhoee, rhoe, rhop, rhow, rhoww,
2792: & rhouee, rhoue, rhoup, rhouw, rhouww,
2793: & ergee, erge, ergp, ergw, ergww,
2794: & jerg,xold)
2795: c
2796: c This function computes the residual
2797: c for the 1-D energy equation
2798: c
2799: c
2800: implicit none
2802: include 'comd.h'
2803: include 'tube.h'
2804: c
2805: c input variables
2806: c
2807: double precision rhoee, rhoe, rhop, rhow, rhoww
2808: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2809: double precision ergee, erge, ergp, ergw, ergww
2810: double precision xold(mx*neq)
2811: c
2812: integer jerg
2813: c
2814: c local variables
2815: c
2816: double precision theta1
2817: integer je
2818: c
2819: c new
2820: c
2821: double precision velfw, velfe
2822: double precision vele,velp,velw
2823: double precision fluxe, fluxw
2824: double precision uepe, uepw
2825: double precision ue, up, uw
2826: double precision see, sep, sew
2827: double precision seef, sewf
2828: double precision upe, upw
2829: double precision presse, pressw
2830: double precision source
2831: double precision te, tp, tw
2832: double precision tef, twf, ref, rwf
2833: double precision kef, kwf
2834: double precision hflxe, hflxw
2835: c
2836: c old
2837: c
2838: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2839: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2840: double precision ergoee, ergoe, ergop, ergow, ergoww
2841: double precision velfow, velfoe
2842: double precision veloe,velop,velow
2843: double precision fluxoe, fluxow
2844: double precision uepoe, uepow
2845: double precision uoe, uop, uow
2846: double precision seoe, seop, seow
2847: double precision seoef, seowf
2848: double precision upoe, upow
2849: double precision pressoe, pressow
2850: double precision sourceo
2851: double precision toe, top, tow
2852: double precision toef, towf, roef, rowf
2853: double precision koef, kowf
2854: double precision hflxoe, hflxow
2855: c
2856: c functions
2857: c
2858: double precision godunov2, eos
2859: double precision upwind, fluxlim
2861: c
2862: c
2863: c ******************************************************************
2864: c
2865: c
2866: c
2867: je = (neq*jerg)
2869: c########################
2870: c
2871: c NEW
2872: c
2873: c########################
2875: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
2876: & rhouee, rhoue, rhoup, rhouw, rhouww,
2877: & ergee, erge, ergp, ergw, ergww,
2878: & vele, velp, velw, jerg)
2880: pressw = eos(rhow, rhouw, ergw)
2881: presse = eos(rhoe, rhoue, erge)
2883: uw = rhouw / rhow
2884: up = rhoup / rhop
2885: ue = rhoue / rhoe
2887: upw = uw * pressw
2888: upe = ue * presse
2890: velfe = 0.5d+0 * (vele + velp)
2891: velfw = 0.5d+0 * (velw + velp)
2893: if (ihod .eq. 1) then
2895: uepe = upwind(ergp,erge,velfe)
2896: uepw = upwind(ergw,ergp,velfw)
2898: elseif (ihod .eq. 2) then
2900: uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2901: uepw = fluxlim(ergww,ergw,ergp,erge,velfw)
2903: endif
2905: if (ihod .eq. 3) then
2906: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
2907: & rhouw,rhoup,rhoue,rhouee,
2908: & ergw, ergp, erge, ergee,3)
2909: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
2910: & rhouww,rhouw,rhoup,rhoue,
2911: & ergww, ergw, ergp, erge,3)
2912: else
2913: fluxe = (dt/dx) * ( uepe + (0.5d+0*upe) )
2914: fluxw = (dt/dx) * ( uepw + (0.5d+0*upw) )
2915: endif
2916: c
2917: c radiation
2918: c
2919: if (kappa0 .eq. 0.0d+0) then
2920: source = 0.0d+0
2921: else
2923: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2924: sep = (ergp/rhop) - (0.5d+0 * up * up)
2925: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2927: seef = 0.5d+0 * (see + sep)
2928: sewf = 0.5d+0 * (sew + sep)
2930: te = see / csubv
2931: tp = sep / csubv
2932: tw = sew / csubv
2934: tef = seef / csubv
2935: twf = sewf / csubv
2937: ref = 0.5d+0 * (rhoe + rhop)
2938: rwf = 0.5d+0 * (rhow + rhop)
2940: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2941: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2943: if (wilson) then
2944: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2945: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2946: endif
2948: if ( debug .and. (kef .gt. 1.0d+10) ) then
2949: write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2950: endif
2951: if ( debug .and. (kwf .gt. 1.0d+10) ) then
2952: write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2953: endif
2955: hflxe = kef * (te - tp) / dx
2956: hflxw = kwf * (tp - tw) / dx
2958: source = (dt/dx) * (hflxe - hflxw)
2960: endif
2962: c########################
2963: c
2964: c OLD
2965: c
2966: c########################
2968: call Setpbc(jerg,xold,
2969: & rhooee, rhooe, rhoop, rhoow, rhooww,
2970: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2971: & ergoee, ergoe, ergop, ergow, ergoww,
2972: & veloe, velop, velow)
2974: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
2975: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
2976: & ergoee, ergoe, ergop, ergow, ergoww,
2977: & veloe, velop, velow, jerg)
2979: pressow = eos(rhoow, rhouow, ergow)
2980: pressoe = eos(rhooe, rhouoe, ergoe)
2982: uow = rhouow / rhoow
2983: uop = rhouop / rhoop
2984: uoe = rhouoe / rhooe
2986: upow = uow * pressow
2987: upoe = uoe * pressoe
2989: velfoe = 0.5d+0 * (veloe + velop)
2990: velfow = 0.5d+0 * (velow + velop)
2993: if (ihod .eq. 1) then
2995: uepoe = upwind(ergop,ergoe,velfoe)
2996: uepow = upwind(ergow,ergop,velfow)
2998: elseif (ihod .eq. 2) then
3000: uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
3001: uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)
3003: endif
3005: if (ihod .eq. 3) then
3006: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3007: & rhouow,rhouop,rhouoe,rhouoee,
3008: & ergow, ergop, ergoe, ergoee,3)
3009: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3010: & rhouoww,rhouow,rhouop,rhouoe,
3011: & ergoww, ergow, ergop, ergoe,3)
3012: else
3013: fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
3014: fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
3015: endif
3017: c
3018: c old radiation
3019: c
3020: if (kappa0 .eq. 0.0d+0) then
3021: sourceo = 0.0d+0
3022: else
3024: seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
3025: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
3026: seow = (ergow/rhoow) - (0.5d+0 * uow * uow)
3028: seoef = 0.5d+0 * (seoe + seop)
3029: seowf = 0.5d+0 * (seow + seop)
3031: toe = seoe / csubv
3032: top = seop / csubv
3033: tow = seow / csubv
3035: toef = seoef / csubv
3036: towf = seowf / csubv
3038: roef = 0.5d+0 * (rhooe + rhoop)
3039: rowf = 0.5d+0 * (rhoow + rhoop)
3041: koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
3042: kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)
3044: if (wilson) then
3045: koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
3046: kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
3047: endif
3049: if ( debug .and. (koef .gt. 1.0d+10) ) then
3050: write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
3051: endif
3052: if ( debug .and. (kowf .gt. 1.0d+10) ) then
3053: write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
3054: endif
3056: hflxoe = koef * (toe - top) / dx
3057: hflxow = kowf * (top - tow) / dx
3059: sourceo = (dt/dx) * (hflxoe - hflxow)
3061: endif
3064: c########################
3065: c
3066: c FUNCTION
3067: c
3068: c########################
3070: CVAM
3071: if (probnum .eq. 3) then
3072: fluxe = 0.0d+0
3073: fluxw = 0.0d+0
3074: fluxoe = 0.0d+0
3075: fluxow = 0.0d+0
3076: endif
3077: CVAM
3079: theta1 = 1.0d+0 - theta
3080: energy = (ergp - xold(je))
3081: & + ( theta * ( (fluxe - fluxw ) - source ) )
3082: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3084: if (debug) then
3085: write(*,*)
3086: write(*,*) 'energy(',jerg,') = ',energy
3087: write(*,*)
3088: write(*,*) fluxe,fluxw
3089: write(*,*) fluxoe,fluxow
3090: write(*,*) source,sourceo
3091: write(*,*)
3092: endif
3094: return
3095: end
3096: double precision function eos(r,ru,e)
3098: implicit none
3100: double precision r,ru,e
3102: double precision se, u
3104: integer ierr
3106: logical debug
3108: include "tube.h"
3110: debug = .false.
3112: if (debug) then
3113: write(*,*)
3114: write(*,*) 'in eos r,ru,e'
3115: write(*,*) r,ru,e
3116: write(*,*)
3117: endif
3119: u = ru/r
3121: se = (e/r) - (0.5d+0 * u * u)
3122: eos = (gamma - 1.0d+0) * r * se
3124: if (eos .lt. 0.0d+0) then
3125: write(*,*)
3126: write(*,*) 'eos = ',eos
3127: write(*,*) 'gamma = ',gamma
3128: write(*,*) 'r = ',r
3129: write(*,*) 'se = ',se
3130: write(*,*) 'e = ',e
3131: write(*,*) 'u = ',u
3132: write(*,*) 'ru = ',ru
3133: call PetscFinalize(ierr)
3134: write(*,*)
3135: stop
3136: endif
3138: if (debug) then
3139: write(*,*)
3140: write(*,*) 'in eos u,se,eos'
3141: write(*,*) u,se,eos
3142: write(*,*)
3143: endif
3146: return
3147: end
3148: subroutine eval2
3150: implicit none
3152: double precision prat, grat, xnum, xdenom
3155: logical debug
3157: include 'tube.h'
3159: debug = .false.
3161: prat = p2/p1
3162: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
3164: xnum = grat + prat
3165: xdenom = 1.0d+0 + (prat * grat)
3166:
3167: e2 = e1 * prat * (xnum/xdenom)
3168:
3171: if (debug) then
3172: write(*,*)
3173: write(*,*) 'e1 = ',e1
3174: write(*,*) 'e2 = ',e2
3175: endif
3177: return
3178: end
3179: subroutine exact0
3181: implicit none
3183: double precision tol, xn
3184: double precision shockp, fprime
3186: integer maxnewt, niter
3188: logical found, debug
3190: include 'tube.h'
3192: debug = .false.
3194: tol = 1.0d-10
3196: maxnewt = 40
3197:
3198: a1 = sqrt(gamma*p1/r1)
3199: a4 = sqrt(gamma*p4/r4)
3203: found = .false.
3204: niter = 0
3206: xn = 0.5d+0 * (p1 + p4)
3208: 10 if ( (.not. found) .and. (niter .le. maxnewt) ) then
3210: niter = niter + 1
3212: xn = xn - (shockp(xn) / fprime(xn))
3214: if (debug) then
3215: write(*,*) niter,shockp(xn),xn
3216: endif
3218: if ( abs(shockp(xn)) .lt. tol ) then
3219: found = .true.
3220: endif
3222: goto 10
3224: endif
3226: if (.not. found) then
3228: write(*,*) 'newton failed'
3229: write(*,*) xn,shockp(xn)
3230: stop
3232: endif
3234: p2 = xn
3237: if (debug) then
3238: write(*,*)
3239: write(*,*) 'p1 = ',p1
3240: write(*,*) 'p2 = ',p2
3241: write(*,*) 'p4 = ',p4
3242: write(*,*)
3243: endif
3245: return
3246: end
3247: double precision function flux(r,ru,e,eqn)
3248: c23456789012345678901234567890123456789012345678901234567890123456789012
3249: c
3250: c function flux
3251: c
3252: c This function computes the flux at a face
3253: c
3254: c23456789012345678901234567890123456789012345678901234567890123456789012
3257: c#######################################################################
3259: implicit none
3261: include 'comd.h'
3262: include 'tube.h'
3264: double precision r, ru, e
3266: integer eqn
3268: double precision p,u
3270: double precision eos
3273: c#######################################################################
3275: p = eos(r,ru,e)
3276: u = ru/r
3278: if (eqn .eq. 1) then
3279: flux = ru
3280: elseif (eqn .eq. 2) then
3281: flux = (u * ru) + p
3282: else
3283: flux = u * (e + p)
3284: endif
3286: return
3287: end
3288: double precision function fluxlim(fww,fw,fe,fee,vp)
3289: c23456789012345678901234567890123456789012345678901234567890123456789012
3290: c
3291: c function fluxlim
3292: c
3293: c this function computes the flux limited quick face value
3294: c
3295: c23456789012345678901234567890123456789012345678901234567890123456789012
3298: c#######################################################################
3300: implicit none
3302: double precision fww, fw, fe, fee, vp
3304: double precision fd, fc, fu
3306: double precision f1, f2, f3, f4, fhod, beta, flc
3308: double precision med, quick
3309:
3310: logical limit
3312: c#######################################################################
3314: limit = .true.
3316: if (vp .gt. 0.0d+0) then
3317: fd = fe
3318: fc = fw
3319: fu = fww
3320: else
3321: fd = fw
3322: fc = fe
3323: fu = fee
3324: endif
3326: fhod = quick(fd,fc,fu)
3328: if (limit) then
3330: beta = 0.25d+0
3331: flc = 4.0d+0
3333: f1 = fc
3334: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3335: f3 = fu + ( flc * (fc - fu) )
3336: f4 = med(f1,f2,f3)
3337: fluxlim = vp * med(f1,f4,fhod)
3339: else
3341: fluxlim = vp * fhod
3343: endif
3345: return
3346: end
3347: double precision function fluxlim2(fww,fw,fe,fee,vp)
3348: c23456789012345678901234567890123456789012345678901234567890123456789012
3349: c
3350: c function fluxlim2
3351: c
3352: c this function computes the flux limited quick face value
3353: c
3354: c23456789012345678901234567890123456789012345678901234567890123456789012
3357: c#######################################################################
3359: implicit none
3361: double precision fww, fw, fe, fee, vp
3363: double precision fd, fc, fu
3365: double precision f1, f2, f3, f4, fhod, beta, flc
3367: double precision med, quick
3368:
3369: logical limit, debug
3371: c#######################################################################
3373: debug = .false.
3375: if (debug) then
3376: write(*,*)
3377: write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3378: write(*,*) fee,fe,fw,fww
3379: write(*,*)
3380: endif
3382: limit = .true.
3384: if (vp .gt. 0.0d+0) then
3385: fd = fe
3386: fc = fw
3387: fu = fww
3388: else
3389: fd = fw
3390: fc = fe
3391: fu = fee
3392: endif
3394: fhod = quick(fd,fc,fu)
3396: if (limit) then
3398: beta = 0.25d+0
3399: flc = 4.0d+0
3401: f1 = fc
3402: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3403: f3 = fu + ( flc * (fc - fu) )
3404: f4 = med(f1,f2,f3)
3405: fluxlim2 = med(f1,f4,fhod)
3407: else
3409: fluxlim2 = fhod
3411: endif
3413: return
3414: end
3415: double precision function fprime(x)
3417: implicit none
3419: double precision x, eps
3420: double precision shockp
3422: eps = 1.0d-8
3424: fprime = ( shockp(x+eps) - shockp(x) ) / eps
3426: return
3427: end
3428: double precision function godent(uhl, uhr, hl, hr, eqn)
3429: c23456789012345678901234567890123456789012345678901234567890123456789012
3430: c
3431: c function godent
3432: c
3433: c this function computes the roe/godunov face value plus entropy fix
3434: c
3435: c23456789012345678901234567890123456789012345678901234567890123456789012
3438: c#######################################################################
3440: implicit none
3442: include 'comd.h'
3444: double precision uhl, uhr, hl, hr, ht, uht
3445: double precision lamr1, lamr2, laml1, laml2
3446: double precision deltal1, deltal2
3448: integer eqn
3450: double precision sum
3452: double precision flux
3454: integer i, j
3456:
3458: c#######################################################################
3460: if (debug) then
3461: write(*,*) 'in godent eqn = ',eqn
3462: endif
3464: c do i = 1,neq
3465: c fr(i) = flux(uhr,hr,i)
3466: c fl(i) = flux(uhl,hl,i)
3467: c enddo
3469: deltau(1) = hr - hl
3470: deltau(2) = uhr - uhl
3472: call roestat(uhl, uhr, hl,hr,ht, uht)
3474: call eigen(ht,uht)
3476: do i = 1,neq
3477: sum = 0.0d+0
3478: do j = 1,neq
3479: sum = sum + ( rinv(i,j) * deltau(j) )
3480: enddo
3481: alpha(i) = sum
3482: enddo
3484: deltal1 = 0.0d+0
3485: deltal2 = 0.0d+0
3487: c call eigene(hr,uhr,lamr1, lamr2)
3488: c call eigene(hl,uhl,laml1, laml2)
3489: c
3490: c 1st eigen
3491: c
3492: if ( (laml1 .lt. 0.0d+0) .and.
3493: & (lamr1 .gt. 0.0d+0) ) then
3495: CVAM deltal1 = 4.0d+0 * (lamr1 - laml1)
3496: deltal1 = 4.0d+0 * (lamr1 - laml1) + 1.0d-2
3498: if ( abs(eigval(1)) .lt. (0.5d+0 * deltal1) ) then
3499: eigval(1) = ( ( (eigval(1) ** 2) / deltal1 )
3500: & + ( 0.25d+0 * deltal1 ) )
3501: endif
3503: endif
3504: c
3505: c 2nd eigen
3506: c
3507: if ( (laml2 .lt. 0.0d+0) .and.
3508: & (lamr2 .gt. 0.0d+0) ) then
3510: deltal2 = 4.0d+0 * (lamr2 - laml2)
3512: if ( abs(eigval(2)) .lt. (0.5d+0 * deltal2) ) then
3513: eigval(2) = ( ( (eigval(2) ** 2) / deltal2 )
3514: & + ( 0.25d+0 * deltal2 ) )
3515: endif
3517: endif
3519: if (debug) then
3520: write(*,*)
3521: write(*,*) 'godent debug'
3522: write(*,*) laml1, laml2, lamr1, lamr2
3523: write(*,*) deltal1, deltal2
3524: write(*,*)
3525: endif
3527: do i = 1,neq
3528: sum = 0.0d+0
3529: do j = 1,neq
3530: sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3531: enddo
3532: xnumdif(i) = sum
3533: enddo
3535: do i = 1,neq
3536: froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3537: enddo
3539: godent = froe(eqn)
3541: return
3542: end
3543: double precision function godunov(uhl, uhr, hl, hr, eqn)
3544: c23456789012345678901234567890123456789012345678901234567890123456789012
3545: c
3546: c function godunov
3547: c
3548: c this function computes the roe/godunov face value
3549: c
3550: c23456789012345678901234567890123456789012345678901234567890123456789012
3553: c#######################################################################
3555: implicit none
3557: include 'comd.h'
3559: double precision uhl, uhr, hl, hr, ht, uht
3561: integer eqn
3563: double precision sum
3565: double precision flux
3567: integer i, j
3568:
3570: c#######################################################################
3572: c do i = 1,neq
3573: c fr(i) = flux(uhr,hr,i)
3574: c fl(i) = flux(uhl,hl,i)
3575: c enddo
3577: deltau(1) = hr - hl
3578: deltau(2) = uhr - uhl
3580: call roestat(uhl, uhr, hl,hr,ht, uht)
3582: call eigen(ht,uht)
3584: do i = 1,neq
3585: sum = 0.0d+0
3586: do j = 1,neq
3587: sum = sum + ( rinv(i,j) * deltau(j) )
3588: enddo
3589: alpha(i) = sum
3590: enddo
3592: do i = 1,neq
3593: sum = 0.0d+0
3594: do j = 1,neq
3595: sum = sum - ( 0.5d+0*alpha(j)*eigvec(i,j)*abs(eigval(j)) )
3596: enddo
3597: xnumdif(i) = sum
3598: enddo
3600: do i = 1,neq
3601: froe(i) = ( 0.5d+0 * (fr(i) + fl(i)) ) + xnumdif(i)
3602: enddo
3604: godunov = froe(eqn)
3607: return
3608: end
3609: double precision function godunov2
3610: & (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3611: c23456789012345678901234567890123456789012345678901234567890123456789012
3612: c
3613: c function godunov2
3614: c
3615: c this function computes the roe/godunov2 face value
3616: c
3617: c23456789012345678901234567890123456789012345678901234567890123456789012
3620: c#######################################################################
3622: implicit none
3624: include 'comd.h'
3625: include 'tube.h'
3627: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3629: integer eqn
3631: double precision rrg, rlg, rurg, rulg, erg, elg
3633: double precision godunov, godent, hlle
3636: c#######################################################################
3638: if (gorder .eq. 1) then
3639: rrg = rr
3640: rlg = rl
3641: rurg = rur
3642: rulg = rul
3643: erg = er
3644: elg = el
3645: else
3646: call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
3647: & rrg, rlg,rurg, rulg, erg, elg)
3648: endif
3650: CVAM if (ientro .eq. 0) then
3651: CVAM godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3652: CVAM elseif(ientro .eq. 1) then
3653: CVAM godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3654: CVAM else
3655: godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3656: CVAM endif
3659: return
3660: end
3661: double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3662: c23456789012345678901234567890123456789012345678901234567890123456789012
3663: c
3664: c function hlle
3665: c
3666: c this function computes the roe/hlle face value
3667: c
3668: c23456789012345678901234567890123456789012345678901234567890123456789012
3671: c#######################################################################
3673: implicit none
3675: include 'comd.h'
3676: include 'tube.h'
3678: double precision rrg,rlg,rurg,rulg,erg,elg
3679: integer eqn
3681: double precision laml1, laml2, laml3
3682: double precision lamr1, lamr2, lamr3
3683: double precision sl, sr
3686: double precision flux
3688: integer i, j, ispeed
3689:
3691: c#######################################################################
3693: ispeed = 1
3695: do i = 1,neq
3696: fr(i) = flux(rrg,rurg,erg,i)
3697: fl(i) = flux(rlg,rulg,elg,i)
3698: enddo
3700: deltau(1) = rrg - rlg
3701: deltau(2) = rurg - rulg
3702: deltau(3) = erg - elg
3704: CVAM call roestat(uhl,uhr,hl,hr,ht,uht)
3706: CVAM call eigene(ht,uht,lamt1, lamt2)
3707: call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3708: call eigene(rlg,rulg,elg,laml1,laml2,laml3)
3710: CVAM if (ispeed .eq. 1) then
3711: CVAM sl = min(laml1,lamt1)
3712: CVAM sr = max(lamt2,lamr2)
3713: CVAM else
3714: sl = min(laml1,lamr1)
3715: sr = max(laml3,lamr3)
3716: CVAM endif
3719: do i = 1,neq
3720: froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )
3721: & / (sr-sl)
3722: enddo
3724: hlle = froe(eqn)
3727: return
3728: end
3729: double precision function med(x1,x2,x3)
3730: c23456789012345678901234567890123456789012345678901234567890123456789012
3731: c
3732: c function med
3733: c
3734: c this function computes the median of three numbers
3735: c
3736: c23456789012345678901234567890123456789012345678901234567890123456789012
3739: c#######################################################################
3741: implicit none
3743: double precision x1, x2, x3
3744: double precision xhi, xlo
3746: c#######################################################################
3748: xhi = max(x1,x2,x3)
3749: xlo = min(x1,x2,x3)
3751: med = x1 + x2 + x3 - xhi - xlo
3753: return
3754: end
3755: double precision function mom
3756: & (rhoee, rhoe, rhop, rhow, rhoww,
3757: & rhouee, rhoue, rhoup, rhouw, rhouww,
3758: & ergee, erge, ergp, ergw, ergww,
3759: & jmom,xold)
3760: c
3761: c This function computes the residual
3762: c for the 1-D momentum equation
3763: c
3764: c
3765: implicit none
3767: include 'comd.h'
3768: include 'tube.h'
3769: c
3770: c input variables
3771: c
3772: double precision rhoee, rhoe, rhop, rhow, rhoww
3773: double precision rhouee, rhoue, rhoup, rhouw, rhouww
3774: double precision ergee, erge, ergp, ergw, ergww
3775: double precision xold(mx*neq)
3776: c
3777: integer jmom
3778: c
3779: c local variables
3780: c
3781: double precision theta1
3782: integer jru
3783: c
3784: c new
3785: c
3786: double precision velfw, velfe
3787: double precision vele,velp,velw
3788: double precision fluxe, fluxw
3789: double precision uurhoe, uurhow
3790: double precision pressee, presse, pressp,pressw, pressww
3791: double precision rupee, rupe, rupp, rupw, rupww
3792: double precision uee, ue, up, uw, uww
3793: double precision source
3794: c
3795: c old
3796: c
3797: double precision rhooee, rhooe, rhoop, rhoow, rhooww
3798: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3799: double precision ergoee, ergoe, ergop, ergow, ergoww
3800: double precision velfow, velfoe
3801: double precision veloe,velop,velow
3802: double precision fluxoe, fluxow
3803: double precision uurhooe, uurhoow
3804: double precision pressoee, pressoe, pressop, pressow, pressoww
3805: double precision rupoee, rupoe, rupop, rupow, rupoww
3806: double precision uoee, uoe, uop, uow, uoww
3807: double precision sourceo
3809: double precision eps
3810: c
3811: c functions
3812: c
3813: double precision godunov2, eos
3814: double precision upwind, fluxlim
3815: c
3816: c
3817: c ******************************************************************
3818: c
3819: c
3820: eps = 1.0d-32
3821: c
3822: jru = (neq*jmom) - 1
3824: c########################
3825: c
3826: c NEW
3827: c
3828: c########################
3830: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,
3831: & rhouee, rhoue, rhoup, rhouw, rhouww,
3832: & ergee, erge, ergp, ergw, ergww,
3833: & vele, velp, velw, jmom)
3835: presse = eos(rhoe, rhoue, erge )
3836: pressw = eos(rhow, rhouw, ergw )
3838: velfe = 0.5d+0 * (vele + velp)
3839: velfw = 0.5d+0 * (velw + velp)
3841: if (ihod .eq. 1) then
3843: uurhoe = upwind(rhoup,rhoue,velfe)
3844: uurhow = upwind(rhouw,rhoup,velfw)
3846: elseif (ihod .eq. 2) then
3848: uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3849: uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)
3851: endif
3853: if (ihod .eq. 3) then
3854: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,
3855: & rhouw,rhoup,rhoue,rhouee,
3856: & ergw, ergp, erge, ergee,2)
3857: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe,
3858: & rhouww,rhouw,rhoup,rhoue,
3859: & ergww, ergw, ergp, erge,2)
3860: else
3861: fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3862: fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3863: endif
3866: source = 0.0d+0
3868: c########################
3869: c
3870: c OLD
3871: c
3872: c########################
3874: call Setpbc(jmom,xold,
3875: & rhooee, rhooe, rhoop, rhoow, rhooww,
3876: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
3877: & ergoee, ergoe, ergop, ergow, ergoww,
3878: & veloe, velop, velow)
3880: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
3881: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
3882: & ergoee, ergoe, ergop, ergow, ergoww,
3883: & veloe, velop, velow, jmom)
3885: pressoe = eos(rhooe, rhouoe, ergoe)
3886: pressow = eos(rhoow, rhouow, ergow)
3888: velfoe = 0.5d+0 * (veloe + velop)
3889: velfow = 0.5d+0 * (velow + velop)
3891: if (ihod .eq. 1) then
3893: uurhooe = upwind(rhouop,rhouoe,velfoe)
3894: uurhoow = upwind(rhouow,rhouop,velfow)
3896: elseif (ihod .eq. 2) then
3898: uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3899: uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)
3901: endif
3903: if (ihod .eq. 3) then
3904: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,
3905: & rhouow,rhouop,rhouoe,rhouoee,
3906: & ergow, ergop, ergoe, ergoee,2)
3907: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,
3908: & rhouoww,rhouow,rhouop,rhouoe,
3909: & ergoww, ergow, ergop, ergoe,2)
3910: else
3911: fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3912: fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3913: endif
3915: sourceo = 0.0d+0
3918: c########################
3919: c
3920: c FUNCTION
3921: c
3922: c########################
3924: theta1 = 1.0d+0 - theta
3925: mom = (rhoup - xold(jru))
3926: & + ( theta * ( (fluxe - fluxw ) - source ) )
3927: & + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3928: CVAM
3929: if (probnum .eq. 3) then
3930: mom = 0.0d+0
3931: endif
3932: CVAM
3933: if (debug) then
3934: write(*,*)
3935: write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3936: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3937: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3938: write(*,*) 'presse = ',presse,'pressw = ',pressw
3939: write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3940: write(*,*)
3941: endif
3943: return
3944: end
3945: double precision function quick(fd, fc, fu)
3946: c23456789012345678901234567890123456789012345678901234567890123456789012
3947: c
3948: c function quick
3949: c
3950: c this function computes the quick face value
3951: c
3952: c23456789012345678901234567890123456789012345678901234567890123456789012
3955: c#######################################################################
3957: implicit none
3959: double precision fd, fc, fu
3961: c#######################################################################
3963: quick = ( (3.0d+0 * fd) + (6.0d+0 * fc) - fu ) / 8.0d+0
3965: return
3966: end
3967: double precision function rexact(x,t)
3969: implicit none
3971: double precision x,t
3972: double precision xot, head, tail, contact, ufan
3973: double precision xpow, grat, urat
3974: double precision uexact
3977: logical debug
3979: include 'tube.h'
3981: debug = .false.
3984: if (t .le. 0.0d+0) then
3985: if (x .gt. 0.0d+0) then
3986: rexact = r1
3987: else
3988: rexact = r4
3989: endif
3990: else
3992: xot = x/t
3993: head = -a4
3994: tail = v3 - a3
3995: contact = v2
3997: if (xot .lt. head) then
3998: rexact = r4
3999: elseif (xot .gt. sspd) then
4000: rexact = r1
4001: elseif (xot .gt. contact) then
4002: rexact = r2
4003: elseif (xot .gt. tail) then
4004: rexact = r3
4005: else
4006: ufan = uexact(x,t)
4007: grat = (gamma - 1.0d+0) / 2.0d+0
4008: xpow = 1.0d+0 / grat
4009: urat = ufan / a4
4010: rexact = r4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
4011: endif
4013: endif
4016: if (debug) then
4017: write(*,*)
4018: write(*,*) 'rexact(',x,',',t,') = ',rexact
4019: write(*,*)
4020: endif
4022: return
4023: end
4024: subroutine roestat(uhl, uhr, hl,hr,ht,uht)
4025: c23456789012345678901234567890123456789012345678901234567890123456789012
4026: c
4027: c subroutine roestat
4028: c
4029: c This subroutine computes the roe state at a face
4030: c
4031: c23456789012345678901234567890123456789012345678901234567890123456789012
4034: c#######################################################################
4036: implicit none
4038: include 'comd.h'
4040: double precision uhl, uhr, hl, hr, ht, uht
4042: double precision ul, ur, shl, shr, xnum, xdenom
4043:
4046: c#######################################################################
4048: ul = uhl / hl
4049: ur = uhr / hr
4051: shl = sqrt(hl)
4052: shr = sqrt(hr)
4054: xnum = (shl * ul) + (shr * ur)
4055: xdenom = shl + shr
4057: ht = 0.5d+0 * (hl + hr)
4058: uht = ht * ( xnum / xdenom )
4060: return
4061: end
4062: subroutine rval2
4064: implicit none
4066: double precision prat, grat, xnum, xdenom
4069: logical debug
4071: include 'tube.h'
4073: debug = .false.
4075: prat = p2/p1
4076: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
4078: xnum = 1.0d+0 + (grat * prat)
4079: xdenom = grat + prat
4080:
4081: r2 = r1 * (xnum/xdenom)
4082:
4085: if (debug) then
4086: write(*,*)
4087: write(*,*) 'r1 = ',r1
4088: write(*,*) 'r2 = ',r2
4089: endif
4091: return
4092: end
4093: subroutine secondq
4094: & (rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,
4095: & rrg, rlg,rurg, rulg, erg, elg)
4096: c23456789012345678901234567890123456789012345678901234567890123456789012
4097: c
4098: c subroutine secondq
4099: c
4100: c this subroutine computes the second order (based on quick) left
4101: c and right states for the godunov solver.
4102: c
4103: c23456789012345678901234567890123456789012345678901234567890123456789012
4106: c#######################################################################
4108: implicit none
4110: include 'comd.h'
4112: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
4113: double precision rrg, rlg,rurg, rulg, erg, elg
4117: double precision veld, ull,ul,ur,urr, ulg, urg
4119: double precision fluxlim2
4122: c#######################################################################
4124: c
4125: c compute the velocities
4126: c
4127: ull = rull/rll
4128: ul = rul /rl
4129: ur = rur /rr
4130: urr = rurr/rrr
4132: c
4133: c compute the left state first
4134: c
4135: veld = 1.0d+0
4137: rlg = fluxlim2(rll,rl,rr,rrr,veld)
4138: ulg = fluxlim2(ull,ul,ur,urr,veld)
4139: rulg = rlg * ulg
4140: elg = fluxlim2(ell,el,er,err,veld)
4141: c
4142: c now compute the right state
4143: c
4144: veld = -1.0d+0
4146: rrg = fluxlim2(rll,rl,rr,rrr,veld)
4147: urg = fluxlim2(ull,ul,ur,urr,veld)
4148: rurg = rrg * urg
4149: erg = fluxlim2(ell,el,er,err,veld)
4153: return
4154: end
4155: double precision function shockp(x)
4157: implicit none
4159: double precision x
4160: double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
4161: logical debug
4163: include 'tube.h'
4165: debug = .false.
4168: if (debug) then
4169: write(*,*)
4170: write(*,*) 'gamma = ',gamma
4171: write(*,*) 'a1 = ',a1
4172: write(*,*) 'a4 = ',a4
4173: write(*,*) 'p1 = ',p1
4174: write(*,*) 'p2 = ',x
4175: write(*,*)
4176: endif
4178: xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
4179: xdenom = sqrt ( 2.0d+0 * gamma * ( (2.0d+0*gamma)
4180: & + (gamma + 1.0d+0) * ((x/p1) - 1) ) )
4181: xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)
4183: shockp = (x/p1)*((1.0d+0-(xnum/xdenom))**xpow) - (p4/p1)
4186: if (debug) then
4187: write(*,*)
4188: write(*,*) 'xnum = ',xnum
4189: write(*,*) 'gamma = ',gamma
4190: write(*,*) 'a1 = ',a1
4191: write(*,*) 'a4 = ',a4
4192: write(*,*) 'p1 = ',p1
4193: write(*,*) 'xdenom = ',xdenom
4194: write(*,*) 'xpow = ',xpow
4195: write(*,*) 'shockp = ',shockp
4196: write(*,*) 'p2 = ',x
4197: write(*,*)
4198: endif
4200: return
4201: end
4202: double precision function uexact(x,t)
4204: implicit none
4206: double precision x,t
4207: double precision xot, head, tail
4210: logical debug
4212: include 'tube.h'
4214: debug = .false.
4216: if (debug) then
4217: write(*,*)
4218: write(*,*) 't = ',t
4219: write(*,*) 'x = ',x
4220: write(*,*) 'a4 = ',a4
4221: write(*,*) 'v3 = ',v3
4222: write(*,*) 'a3 = ',a3
4223: write(*,*)
4224: endif
4226: if (t .le. 0.0d+0) then
4227: uexact = 0.0d+0
4228: else
4230: xot = x/t
4231: head = -a4
4232: tail = v3 - a3
4234: if (xot .lt. head) then
4235: uexact = 0.0d+0
4236: elseif (xot .gt. sspd) then
4237: uexact = 0.0d+0
4238: elseif (xot .gt. tail) then
4239: uexact = v2
4240: else
4241: uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
4242: endif
4244: endif
4247: if (debug) then
4248: write(*,*)
4249: CVAM write(*,*) 'x = ',x,' t = ',t
4250: write(*,*) 'uexact = ',uexact
4251: write(*,*)
4252: endif
4254: return
4255: end
4256: double precision function upwind(fw, fe, vp)
4257: c23456789012345678901234567890123456789012345678901234567890123456789012
4258: c
4259: c function upwind
4260: c
4261: c this function computes the upwind face value
4262: c
4263: c23456789012345678901234567890123456789012345678901234567890123456789012
4266: c#######################################################################
4268: implicit none
4270: double precision fw, fe, vp
4272: c#######################################################################
4274: if (vp .gt. 0.0) then
4275: upwind = vp * fw
4276: else
4277: upwind = vp * fe
4278: endif
4280: return
4281: end
4282: subroutine uval2
4284: implicit none
4286: double precision prat, grat1, grat2, arat, xnum
4289: logical debug
4291: include 'tube.h'
4293: debug = .false.
4295: prat = p2/p1
4296: grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
4297: grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
4298: arat = a1/gamma
4300: xnum = sqrt ( grat2 / (prat + grat1) )
4302: v2 = arat * (prat - 1.0d+0) * xnum
4304: if (debug) then
4305: write(*,*)
4306: write(*,*) 'v2 = ',v2
4307: endif
4309: return
4310: end
4311: subroutine val3
4313: implicit none
4315: double precision prat, rpow, epow, p3t
4318: logical debug
4320: include 'tube.h'
4322: debug = .false.
4325: p3 = p2
4327: prat = p3/p4
4329: rpow = 1.0d+0 / gamma
4331: r3 = r4 * ( prat ** rpow )
4333: epow = (gamma - 1.0d+0) / gamma
4335: e3 = e4 * ( (p3/p4) ** epow )
4337: p3t = (gamma - 1.0d+0) * r3 * e3
4339: a3 = sqrt(gamma*p3/r3)
4341: if (debug) then
4342: write(*,*)
4343: write(*,*) 'a3 = ',a3
4344: write(*,*) 'r3 = ',r3
4345: write(*,*) 'e3 = ',e3
4346: write(*,*) 'p3 = ',p3
4347: write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
4348: write(*,*)
4349: endif
4351: return
4352: end
4353: subroutine wval
4355: implicit none
4357: double precision prat, grat, xnum
4360: logical debug
4362: include 'tube.h'
4364: debug = .false.
4366: prat = p2/p1
4367: grat = (gamma + 1.0d+0) / (2.0d+0 * gamma)
4369: xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
4370:
4371: sspd = a1 * sqrt(xnum)
4372:
4375: if (debug) then
4376: write(*,*)
4377: write(*,*) 'sspd = ',sspd
4378: endif
4380: return
4381: end