Actual source code: ex74f.F90
petsc-3.5.4 2015-05-23
1: program radhyd
2: ! "$Id: ex4f.F,v 1.39 1999/03/10 19:29:25 Vince Mousseau $";
3: !
4: ! Description: This example solves a nonlinear system on 1 processor with SNES.
5: ! We solve the Euler equations in one dimension.
6: ! The command line options include:
7: ! -dt <dt>, where <dt> indicates time step
8: ! -mx <xg>, where <xg> = number of grid points in the x-direction
9: ! -nstep <nstep>, where <nstep> = number of time steps
10: ! -debug <ndb>, where <ndb> = 0) no debug 1) debug
11: ! -pcnew <npc>, where <npc> = 0) no preconditioning 1) rad preconditioning
12: ! -probnum <probnum>, where <probnum> = 1) cyclic Riesner 2) dam break
13: ! -ihod <ihod>, where <ihod> = 1) upwind 2) quick 3) godunov
14: ! -ientro <ientro>, where <ientro> = 0) basic 1) entropy fix 2) hlle
15: ! -theta <theta>, where <theta> = 0-1 0-explicit 1-implicit
16: ! -hnaught <hnaught>, where <hnaught> = height of left side
17: ! -hlo <hlo>, where <hlo> = hieght of right side
18: ! -ngraph <ngraph>, where <ngraph> = number of time steps between graphics
19: ! -damfac <damfac>, where <damfac> = fractional downward change in hight
20: ! -dampit <ndamp>, where <ndamp> = 1 turn damping on
21: ! -gorder <gorder>, where <gorder> = spatial oerder of godunov
22: !
23: !
24: ! --------------------------------------------------------------------------
25: !
26: ! Shock tube example
27: !
28: ! In this example the application context is a Fortran integer array:
29: ! ctx(1) = shell preconditioner pressure matrix contex
30: ! ctx(2) = semi implicit pressure matrix
31: ! ctx(4) = xold - old time values need for time advancement
32: ! ctx(5) = mx - number of control volumes
33: ! ctx(6) = N - total number of unknowns
34: !
35: ! --------------------------------------------------------------------------
37: implicit none
39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: ! Include files
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: !
43: #include finclude/petscdef.h
44: #include finclude/petsc.h
46: #include "ex74fcomd.h"
47: #include "ex74ftube.h"
49: !
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Variable declarations
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: !
54: ! Variables:
55: ! snes - nonlinear solver
56: ! x,r - solution, residual vectors
57: ! J - Jacobian matrix
58: ! its - iterations for convergence
59: ! dt - time step size
60: ! draw - drawing context
61: !
62: PetscFortranAddr ctx(6)
63: integer rank,size,nx,ny
64: SNES snes
65: KSP ksp
66: PC pc
67: Vec x,r
68: PetscViewer view0, view1, view2, view3, view4
69: Mat Psemi
70: integer flg, N, ierr, ngraph
71: integer nstep, ndt, i
72: integer its, lits, totits, totlits
73: integer ndb, npc, ndamp, nwilson, ndtcon
74: double precision plotim
76: double precision krtol,katol,kdtol
77: double precision natol,nrtol,nstol
78: integer kmit,nmit,nmf
81: ! Note: Any user-defined Fortran routines (such as FormJacobian)
82: ! MUST be declared as external.
84: external FormFunction, FormInitialGuess,FormDt,PCRadSetUp, PCRadApply, FormGraph,FormDampit
85: double precision eos
87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: ! Initialize program
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: open (unit=87,file='Dt.out',status='unknown')
94: ! start PETSc
95: !
96: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
97: call PetscOptionsSetValue('-snes_mf','true',ierr)
98: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
99: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
101: if (size .ne. 1) then
102: if (rank .eq. 0) then
103: write(6,*) 'This is a uniprocessor example only!'
104: endif
105: SETERRQ(PETSC_COMM_WORLD,1,0,' ')
106: endif
108: ! Initialize problem parameters
110: debug = .false.
111: dampit = .false.
112: wilson = .true.
113: dtcon = .true.
114: pcnew = .true.
115: dtmax = 1.0d+2
116: dtmin = 1.00d-12
117: dt = 1.0d-2
118: mx = 100
119: nstep = 50
120: probnum = 1
121: gorder = 1
123: tfinal = 1.0d+0
124: tplot = 0.2d+0
125: dtgrow = 1.05d+0
126: tcscal = 0.5d+0
127: hcscal = 0.5d+0
129: ihod = 3
130: ientro = 1
131: theta = 1.00d+0
132: pi = 3.14159d+0
134: zero = 0.0
135: ngraph = 10
137: ndb = 0
138: npc = 1
140: damfac = 0.9d+0
142: gamma = 1.25d+0
143: csubv = 1.0d+0 / (gamma - 1.0d+0)
145: v1 = 0.0d+0
146: v4 = 0.0d+0
148: e1 = 1.0d+0
149: e4 = 1.0d+0
151: r1 = 1.0d+0
152: r4 = 2.0d+0
154: ru1 = r1 * v1
155: ru4 = r4 * v4
157: et1 = r1 * ( (0.5d+0 * v1 * v1) + e1 )
158: et4 = r4 * ( (0.5d+0 * v4 * v4) + e4 )
160: p1 = eos(r1,ru1,et1)
161: p4 = eos(r4,ru4,et4)
163: a1 = sqrt(gamma*p1/r1)
164: a4 = sqrt(gamma*p4/r4)
166: erg0 = 1.0d+2
167: kappa0 = 1.0d+0
168: kappaa = -2.0d+0
169: kappab = 13.0d+0 / 2.0d+0
171: !
172: ! load the command line options
173: !
174: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
175: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
176: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nstep',nstep,flg,ierr)
177: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-debug',ndb,flg,ierr)
178: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-pcnew',npc,flg,ierr)
179: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ihod',ihod,flg,ierr)
180: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ientro',ientro,flg,ierr)
181: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',theta,flg,ierr)
182: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ngraph',ngraph,flg,ierr)
183: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-damfac',damfac,flg,ierr)
184: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dampit',ndamp,flg,ierr)
185: call PetscOptionsGetInt(PETSC_NULL_CHARACTER, '-wilson',nwilson,flg,ierr)
186: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-gorder',gorder,flg,ierr)
187: call PetscOptionsGetInt(PETSC_NULL_CHARACTER, '-probnum',probnum,flg,ierr)
188: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-kappa0',kappa0,flg,ierr)
189: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-erg0',erg0,flg,ierr)
190: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dtcon',ndtcon,flg,ierr)
191: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tfinal',tfinal,flg,ierr)
192: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tplot',tplot,flg,ierr)
193: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dtgrow',dtgrow,flg,ierr)
194: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-tcscal',tcscal,flg,ierr)
195: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-hcscal',hcscal,flg,ierr)
196: call PetscOptionsGetReal(PETSC_NULL_CHARACTER, '-dtmax',dtmax,flg,ierr)
198: if (ndamp .eq. 1) then
199: dampit = .true.
200: endif
202: if (nwilson .eq. 0) then
203: wilson = .false.
204: endif
206: if (ndb .eq. 1) then
207: debug = .true.
208: endif
210: if (npc .eq. 0) then
211: pcnew = .false.
212: endif
214: if (ndtcon .eq. 0) then
215: dtcon = .false.
216: endif
218: !CVAM if (dt .ge. dtmax .or. dt .le. dtmin) then
219: !CVAM if (rank .eq. 0) write(6,*) 'DT is out of range'
220: !CVAM SETERRA(1,0,' ')
221: !CVAM endif
223: N = mx*neq
225: ctx(5) = mx
226: ctx(6) = N
228: if (debug) then
229: write(*,*) 'mx = ',mx
230: endif
234: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: ! Create nonlinear solver context
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
240: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: ! Create vector data structures; set function evaluation routine
242: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244: call MatCreateSeqAIJ(PETSC_COMM_WORLD,mx,mx,10,PETSC_NULL_INTEGER,ctx(2),ierr)
246: if (debug) then
247: call MatGetSize(ctx(2),nx,ny,ierr)
248: write(*,*) 'number of rows = ',nx,' number of col = ',ny
249: endif
250: !
251: ! full size vectors
252: !
253: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,x,ierr)
254: call VecSetFromOptions(x,ierr)
255: call VecDuplicate(x,r,ierr)
256: call VecDuplicate(x,ctx(4),ierr)
257: !
258: ! set grid
259: !
260: dx = 2.0d+0/dfloat(mx)
261: xl0 = -1.0d+0 -(0.5d+0 * dx)
263: if (debug) then
264: write(*,*) 'dx = ',dx
265: endif
266:
268: ! Set function evaluation routine and vector. Whenever the nonlinear
269: ! solver needs to evaluate the nonlinear function, it will call this
270: ! routine.
271: ! - Note that the final routine argument is the user-defined
272: ! context that provides application-specific data for the
273: ! function evaluation routine.
275: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
277: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: ! Customize nonlinear solver; set runtime options
279: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
283: call SNESSetFromOptions(snes,ierr)
284: !
285: ! set the line search function to damp the newton update.
286: !
287: if (dampit) then
288: call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
289: endif
290: !
291: ! get the linear solver info from the nonlinear solver
292: !
294: call SNESGetKSP(snes,ksp,ierr)
295: call KSPGetPC(ksp,pc,ierr)
297: call KSPGetTolerances(ksp,krtol,katol,kdtol,kmit,ierr)
298: call SNESGetTolerances(snes,natol,nrtol,nstol,nmit,nmf,ierr)
300: write(*,*)
301: write(*,*)
302: write(*,*) 'Linear solver'
303: write(*,*)
304: write(*,*) 'rtol = ',krtol
305: write(*,*) 'atol = ',katol
306: write(*,*) 'dtol = ',kdtol
307: write(*,*) 'maxits = ',kmit
308: write(*,*)
309: write(*,*)
310: write(*,*) 'Nonlinear solver'
311: write(*,*)
312: write(*,*) 'rtol = ',nrtol
313: write(*,*) 'atol = ',natol
314: write(*,*) 'stol = ',nstol
315: write(*,*) 'maxits = ',nmit
316: write(*,*) 'max func = ',nmf
317: write(*,*)
318: write(*,*)
320: !
321: ! Build shell based preconditioner if flag set
322: !
323: if (pcnew) then
324: call PCSetType(pc,PCSHELL,ierr)
325: call PCShellSetContext(pc,ctx,ierr)
326: call PCShellSetSetUpCtx(pc,PCRadSetUp,ierr)
327: call PCShellSetApplyCtx(pc,PCRadApply,ctx,ierr)
328: endif
330: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
332: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: ! Evaluate initial guess; then solve nonlinear system.
334: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335: !
336: ! initial counters
337: !
338: time = 0.0d+0
339: plotim = 0.0d+0
340: totits = 0
341: totlits = 0
343: ! Note: The user should initialize the vector, x, with the initial guess
344: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
345: ! to employ an initial guess of zero, the user should explicitly set
346: ! this vector to zero by calling VecSet().
348: call FormInitialGuess(x,ierr)
349: !
350: ! open a window to plot results
351: !
352: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'density',0,0,300,300,view0,ierr)
353: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'velocity',320,0,300,300,view1,ierr)
354: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'total energy',640,0,300,300,view2,ierr)
355: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'temperature',0,360,300,300,view3,ierr)
356: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'pressure',320,360,300,300,view4,ierr)
357: !
358: ! graph initial conditions
359: !
360: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
361: !
362: ! copy x into xold
363: !
364: call VecCopy(x,ctx(4),ierr)
365: call FormDt(snes,x,ctx,ierr)
366: !
367: !################################
368: !
369: ! TIME STEP LOOP BEGIN
370: !
371: !################################
372: !
373: ndt = 0
375: 10 if ( (ndt .le. nstep) .and. ((time + 1.0d-10) .lt. tfinal) ) then
377: if (debug) then
378: write(*,*)
379: write(*,*) 'start of time loop'
380: write(*,*)
381: write(*,*) 'ndt = ',ndt
382: write(*,*) 'nstep = ',nstep
383: write(*,*) 'time = ',time
384: write(*,*) 'tfinal = ',tfinal
385: write(*,*)
386: endif
388: ndt = ndt + 1
389: !
390: ! increment time
391: !
392: time = time + dt
393: plotim = plotim + dt
394: !
395: ! call the nonlinear solver
396: !
397: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
398: !
399: ! get the number of linear iterations used by the nonlinear solver
400: !
401: call SNESGetLinearSolveIterations(snes,lits,ierr)
402: call SNESGetIterationNumber(snes,its,ierr)
404: if (debug) then
405: write(*,*) 'in radhyd ',ndt,'x'
406: call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
407: endif
408: !
409: ! increment the counters
410: !
411: totits = totits + its
412: totlits = totlits + lits
413: !
414: ! Compute new time step
415: !
416: call FormDt(snes,x,ctx,ierr)
417: !
418: ! Draw contour plot of solution
419: !
420: if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
421:
422: plotim = plotim - tplot
425: if (rank .eq. 0) then
426: write(6,100) totits,totlits,ndt,dt,time
427: endif
428: 100 format('Newt = ',i7,' lin =',i7,' step =',i7,' dt = ',e8.3,' time = ',e10.4)
429: !
430: ! graph state conditions
431: !
432: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
434: endif
435: !
436: ! copy x into xold
437: !
438: call VecCopy(x,ctx(4),ierr)
441: goto 10
443: endif
444: !
445: !################################
446: !
447: ! TIME STEP LOOP END
448: !
449: !################################
450: !
452: !
453: ! graph final conditions
454: !
455: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
458: write(*,*)
459: write(*,*)
460: write(*,*) 'total Newton iterations = ',totits
461: write(*,*) 'total linear iterations = ',totlits
462: write(*,*) 'Average Newton per time step = ', dble(totits)/dble(ndt)
463: write(*,*) 'Average Krylov per Newton = ', dble(totlits)/dble(totits)
464: write(*,*)
465: write(*,*)
467: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
468: ! Free work space. All PETSc objects should be destroyed when they
469: ! are no longer needed.
470: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473: call MatDestroy(ctx(2),ierr)
474: call VecDestroy(x,ierr)
475: call VecDestroy(ctx(4),ierr)
476: call VecDestroy(r,ierr)
477: call SNESDestroy(snes,ierr)
478: call PetscViewerDestroy(view0,ierr)
479: call PetscViewerDestroy(view1,ierr)
480: call PetscViewerDestroy(view2,ierr)
481: call PetscViewerDestroy(view3,ierr)
482: call PetscViewerDestroy(view4,ierr)
484: call PCDestroy(ctx(1),ierr)
486: call PetscFinalize(ierr)
488: close(87)
490: stop
491: end
492: subroutine ApplicationDampit(x,deltx,w,ierr)
493: ! ---------------------------------------------------------------------
494: !
495: ! ApplicationDampit - Damps the newton update, called by
496: ! the higher level routine FormDampit().
497: !
498: ! Input Parameter:
499: ! x - current iterate
500: ! deltx - update
501: !
502: ! Output Parameters:
503: ! x - new iterate
504: ! ierr - error code
505: !
506: ! Notes:
507: ! This routine only damps the density. May want to add energy
508: ! in the future
509: !
511: implicit none
513: ! Common blocks:
514: #include "ex74fcomd.h"
516: ! Input/output variables:
517: PetscScalar x(mx*neq), deltx(mx*neq), w(mx*neq)
518: integer ierr
520: ! Local variables:
521: double precision facmin, fac, newx, xmin, se, dse
522: double precision u,en,rn,run
523: integer i, jr, jru, je
525: 0
527: if (debug) then
528: write(*,*) 'begin damping'
529: do i = 1,mx*neq
530: write(*,*)'x(',i,') = ',x(i)
531: enddo
532: write(*,*)
533: do i = 1,mx*neq
534: write(*,*)'deltx(',i,') = ',deltx(i)
535: enddo
536: endif
538: facmin = 1.0d+0
539: !
540: ! set the scale factor
541: !
542: do i=1,mx
543: !
544: ! set pointers
545: !
546: jr = (neq*i) - 2
547: jru = (neq*i) - 1
548: je = (neq*i)
549: !
550: ! make sure dencity stayes positive
551: !
552: newx = x(jr) - deltx(jr)
553: xmin = damfac * x(jr)
555: if (newx .lt. xmin) then
556: fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
557: if (fac .lt. facmin) then
558: if (debug) then
559: write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
560: endif
561: facmin = fac
562: endif
563: endif
564: !
565: ! make sure Total energy stayes positive
566: !
567: newx = x(je) - deltx(je)
568: xmin = damfac * x(je)
570: if (newx .lt. xmin) then
571: fac = (1.0d+0 - damfac)*x(je)/deltx(je)
572: if (fac .lt. facmin) then
573: if (debug) then
574: write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
575: endif
576: facmin = fac
577: endif
578: endif
579: !
580: ! make sure specific internal energy stayes positive
581: !
582:
583: u = x(jru)/x(jr)
584: se = (x(je)/x(jr)) - (0.5d+0 * u * u)
586: en = x(je) - deltx(je)
587: rn = x(jr) - deltx(jr)
588: run = x(jru) - deltx(jru)
590: dse = se - ( (en/rn) - (0.5d+0 * (run/rn) * (run/rn)) )
593: newx = se - dse
594: xmin = damfac * se
596: if (newx .lt. xmin) then
597: fac = (1.0d+0 - damfac) * se / dse
598: if (fac .lt. facmin) then
599: if (debug) then
600: write(*,*) 'se',i, damfac,facmin,fac,se,dse
601: endif
602: facmin = fac
603: endif
604: endif
606: enddo
607: !
608: ! write out warning
609: !
610: if (facmin .lt. 1.0d+0) then
611: write(*,*) 'facmin = ',facmin, damfac,time
612: !
613: ! scale the vector
614: !
615: do i=1,neq*mx
616: w(i) = x(i) - (facmin * deltx(i))
617: enddo
618: else
619: do i=1,neq*mx
620: w(i) = x(i) - deltx(i)
621: enddo
622: endif
624: if (debug) then
625: write(*,*) 'end damping'
626: do i = 1,mx*neq
627: write(*,*) 'w(',i,') = ',w(i)
628: enddo
629: endif
631: return
632: end
633: subroutine ApplicationDt(x,xold,ierr)
634: ! ---------------------------------------------------------------------
635: !
636: ! ApplicationDt - compute CFL numbers. Called by
637: ! the higher level routine FormDt().
638: !
639: ! Input Parameter:
640: ! x - local vector data
641: !
642: ! Output Parameters:
643: ! ierr - error code
644: !
645: ! Notes:
646: ! This routine uses standard Fortran-style computations over a 2-dim array.
647: !
649: implicit none
651: ! Common blocks:
652: #include "ex74fcomd.h"
653: #include "ex74ftube.h"
655: ! Input/output variables:
656: PetscScalar x(mx*neq), xold(mx*neq)
657: integer ierr
659: ! Local variables:
660: integer i, jr, jru, je
661: !
662: ! new
663: !
664: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw
665: double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
666: double precision tcflg, hcflg, dtt, dth
667: double precision te, tp, tw
668: double precision ue, up, uw
669: double precision see, sep, sew
670: !
671: ! old
672: !
673: double precision rhooee, rhooe, rhoop, rhoow, rhooww
674: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
675: double precision ergoee, ergoe, ergop, ergow, ergoww
676: double precision veloe, velop, velow
677: double precision uop, seop, top
678: double precision dtold, dttype
679: !
680: ! functions
681: !
682: double precision eos
684: dtold = dt
686: 0
688: if (debug) then
689: write(*,*) 'in start dt'
690: do i = 1,mx*neq
691: write(*,*)'x(',i,') = ',x(i)
692: enddo
693: write(*,*) 'tfinal = ',tfinal
694: write(*,*) 'time = ',time
695: write(*,*) 'dt = ',dt
696: write(*,*) 'dtmax = ',dtmax
697: endif
699: sndp = -1.0d+20
700: vradn = 0.0d+0
701: vradd = 0.0d+0
703: !
704: !################################
705: !
706: ! loop over all cells begin
707: !
708: !################################
709: !
710: do i=1,mx
711: !
712: ! set pointers
713: !
714: jr = (neq*i) - 2
715: jru = (neq*i) - 1
716: je = (neq*i)
717: !
718: !
719: ! set scalars
720: !
721: call Setpbc(i,x, rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww,vele, velp, velw)
722: !
723: ! compute temperatures
724: !
725: uw = rhouw / rhow
726: up = rhoup / rhop
727: ue = rhoue / rhoe
729: see = (erge/rhoe) - (0.5d+0 * ue * ue)
730: sep = (ergp/rhop) - (0.5d+0 * up * up)
731: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
733: te = see / csubv
734: tp = sep / csubv
735: tw = sew / csubv
736: !
737: ! compute old temperature
738: !
740: call Setpbc(i,xold,rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
742: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww,veloe, velop, velow, i)
744: uop = rhouop / rhoop
746: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
748: top = seop / csubv
749: !
750: ! compute thermal cfl
751: !
752: vradn = vradn + (abs(tp - top)/dt)
753: vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
754: !
755: ! compute hydro cfl
756: !
758: pressp = eos(rhop, rhoup, ergp)
759: sndp = max(sndp,sqrt( (gamma*pressp) / rhop ))
761: enddo
762: !
763: !################################
764: !
765: ! loop over all cells end
766: !
767: !################################
768: !
770: vrad = vradn / vradd
772: tcfl = (vrad * dt) / dx
773: hcfl = (sndp * dt) / dx
775: dtt = max(dx/vrad,1.0d-7)
776: dtt = tcscal * dtt
778: dth = hcscal * dx / sndp
780: if (.not. dtcon) then
781: dt = min (dth,dtt,dt*dtgrow)
782: if (dt .lt. dtmin) then
783: dt = dtmin
784: endif
785: if (dt .gt. dtmax) then
786: dt = dtmax
787: endif
788: if ( (time + dt) .gt. tfinal) then
789: dt = tfinal - time
790: endif
792: if (dt .eq. dth) then
793: dttype = 1.0d+0
794: elseif (dt .eq. dtt) then
795: dttype = 2.0d+0
796: elseif (dt .eq. dtold*dtgrow) then
797: dttype = 3.0d+0
798: elseif (dt .eq. dtmax) then
799: dttype = 4.0d+0
800: elseif (dt .eq. dtmin) then
801: dttype = 5.0d+0
802: elseif (dt .eq. tfinal - time) then
803: dttype = 6.0
804: else
805: dttype = -1.0d+0
806: endif
808: endif
809:
810:
811: write(87,1000) time,dt,dth/hcscal,dtt/tcscal
812: write(88,1000) time,dttype
814: 1000 format(4(2x,e18.9))
816: if (debug) then
817: write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
818: write(*,*) 'dtt = ',dtt,' dth = ',dth
819: write(*,*) 'tfinal = ',tfinal
820: write(*,*) 'time = ',time
821: write(*,*) 'dt = ',dt
822: write(*,*) 'dtmax = ',dtmax
823: write(*,*)
824: endif
826: return
827: end
828: subroutine ApplicationExact(x,ierr)
829: ! ---------------------------------------------------------------------
830: !
831: ! ApplicationExact - Computes exact solution, called by
832: ! the higher level routine FormExact().
833: !
834: ! Input Parameter:
835: ! x - local vector data
836: !
837: ! Output Parameters:
838: ! x - initial conditions
839: ! ierr - error code
840: !
841: ! Notes:
842: ! This routine uses standard Fortran-style computations over a 1-dim array.
843: !
845: implicit none
847: ! Common blocks:
849: #include "ex74fcomd.h"
851: ! Input/output variables:
852: PetscScalar x(mx)
853: integer ierr
855: ! Local variables:
856: integer i
857: double precision xloc
858: PetscScalar rexact
861: ! Set parameters
863: 0
865: do i = 1,mx
867: xloc = xl0 + (dble(i) * dx)
868: x(i) = rexact(xloc,time)
870: enddo
872: return
873: end
874: subroutine ApplicationFunction(x,f,xold,ierr)
875: ! ---------------------------------------------------------------------
876: !
877: ! ApplicationFunction - Computes nonlinear function, called by
878: ! the higher level routine FormFunction().
879: !
880: ! Input Parameter:
881: ! x - local vector data
882: !
883: ! Output Parameters:
884: ! f - local vector data, f(x)
885: ! ierr - error code
886: !
887: ! Notes:
888: ! This routine uses standard Fortran-style computations over a 2-dim array.
889: !
891: implicit none
893: ! Common blocks:
894: #include "ex74fcomd.h"
896: ! Input/output variables:
897: PetscScalar x(mx*neq), f(mx*neq), xold(mx*neq)
898: integer ierr
900: ! Local variables:
901: integer i, jr, jru, je
902: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw
904: double precision cont, energy, mom
906: 0
908: if (debug) then
909: write(*,*) 'in function'
910: do i = 1,mx*neq
911: write(*,*)'x(',i,') = ',x(i)
912: enddo
913: endif
914: !
915: !################################
916: !
917: ! loop over all cells begin
918: !
919: !################################
920: !
921: do i=1,mx
922: !
923: ! set pointers
924: !
925: jr = (neq*i) - 2
926: jru = (neq*i) - 1
927: je = (neq*i)
928: !
929: !
930: ! set scalars
931: !
932: call Setpbc(i,x,rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw)
933: !
934: ! compute functions
935: !
937: f(jr) = cont(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
940: f(jru) = mom(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
943: f(je) = energy(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
945: if (debug) then
946: write(*,*)
947: write(*,*) i,jr,jru,je,'res,r,ru,e'
948: write(*,*) f(jr),f(jru),f(je)
949: write(*,*)
950: endif
952: enddo
953: !
954: !################################
955: !
956: ! loop over all cells end
957: !
958: !################################
959: !
961: if (debug) then
962: write(*,*) 'in function'
963: do i = 1,mx*neq
964: write(*,*) 'f(',i,') = ',f(i)
965: enddo
966: endif
968: return
969: end
970: subroutine ApplicationInitialGuess(x,ierr)
971: ! ---------------------------------------------------------------------
972: !
973: ! ApplicationInitialGuess - Computes initial approximation, called by
974: ! the higher level routine FormInitialGuess().
975: !
976: ! Input Parameter:
977: ! x - local vector data
978: !
979: ! Output Parameters:
980: ! x - initial conditions
981: ! ierr - error code
982: !
983: ! Notes:
984: ! This routine uses standard Fortran-style computations over a 1-dim array.
985: !
987: implicit none
989: ! Common blocks:
991: #include "ex74fcomd.h"
992: #include "ex74ftube.h"
994: ! Input/output variables:
995: PetscScalar x(mx*neq)
996: integer ierr
998: ! Local variables:
999: integer i, j, jr, jru, je
1000: double precision xloc, re, ee, ve
1001: double precision wid, efloor
1002: PetscScalar uexact, rexact, eexact
1005: !VAM efloor = max(1.0d-1,1.0d-3 * erg0)
1006: efloor = 1.0d-1
1007: !VAM wid = max(1.0d-2,dx)
1008: wid = 1.0d-2
1010: ! Set parameters
1012: 0
1014: do i = 1,mx
1016: jr = (neq*i) - 2
1017: jru = (neq*i) - 1
1018: je = (neq*i)
1020: xloc = xl0 + (dble(i) * dx)
1022: if (probnum .eq. 1) then
1023: re = rexact(xloc,zero)
1024: ve = uexact(xloc,zero)
1025: ee = eexact(xloc,zero)
1026: else
1027: re = 1.0d+0
1028: ve = 0.0d+0
1029: ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1030: endif
1032: x(jr) = re
1033: x(jru) = re * ve
1034: x(je) = re * ( (0.5d+0 * ve * ve) + ee )
1036: if (debug) then
1037: write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1038: 100 format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1039: endif
1041: enddo
1043: call exact0
1044: call eval2
1045: call rval2
1046: call wval
1047: call uval2
1048: v3 = v2
1049: call val3
1051: a1 = sqrt(gamma*p1/r1)
1052: a2 = sqrt(gamma*p2/r2)
1053: a3 = sqrt(gamma*p3/r3)
1054: a4 = sqrt(gamma*p4/r4)
1056: write(*,1000) r1,r2,r3,r4
1057: write(*,2000) p1,p2,p3,p4
1058: write(*,3000) e1,e2,e3,e4
1059: write(*,4000) a1,a2,a3,a4
1060: write(*,*)
1062: 1000 format ('rhos ',4(f12.6))
1063: 2000 format ('pressures ',4(f12.6))
1064: 3000 format ('energies ',4(f12.6))
1065: 4000 format ('sound ',4(f12.6))
1068: return
1069: end
1070: subroutine ApplicationXmgr(x,ivar,ierr)
1071: ! ---------------------------------------------------------------------
1072: !
1073: ! ApplicationXmgr - Sets the Xmgr output called from
1074: ! the higher level routine FormXmgr().
1075: !
1076: ! Input Parameter:
1077: ! x - local vector data
1078: !
1079: ! Output Parameters:
1080: ! x - initial conditions
1081: ! ierr - error code
1082: !
1083: ! Notes:
1084: ! This routine uses standard Fortran-style computations over a 1-dim array.
1085: !
1087: implicit none
1089: ! Common blocks:
1091: #include "ex74fcomd.h"
1093: ! Input/output variables:
1094: PetscScalar x(mx)
1095: integer ivar,ierr
1097: ! Local variables:
1098: integer i
1099: double precision xloc, sum
1100: PetscScalar rexact
1101: integer iplotnum(5)
1102: save iplotnum
1103: character*8 grfile
1106: data iplotnum / -1,-1,-1,-1,-1 /
1110: ! Set parameters
1112: iplotnum(ivar) = iplotnum(ivar) + 1
1113: 0
1115: if (ivar .eq. 1) then
1116: write(grfile,4000) iplotnum(ivar)
1117: 4000 format('Xmgrr',i3.3)
1118: elseif (ivar .eq. 2) then
1119: write(grfile,5000) iplotnum(ivar)
1120: 5000 format('Xmgru',i3.3)
1121: elseif (ivar .eq. 3) then
1122: write(grfile,6000) iplotnum(ivar)
1123: 6000 format('Xmgre',i3.3)
1124: elseif (ivar .eq. 4) then
1125: write(grfile,7000) iplotnum(ivar)
1126: 7000 format('Xmgrt',i3.3)
1127: else
1128: write(grfile,8000) iplotnum(ivar)
1129: 8000 format('Xmgrp',i3.3)
1130: endif
1132: open (unit=44,file=grfile,status='unknown')
1134: do i = 1,mx
1136: xloc = xl0 + (dble(i) * dx)
1137: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1138: write(44,1000) xloc, x(i), rexact(xloc,time)
1139: else
1140: write(44,1000) xloc, x(i)
1141: endif
1143: enddo
1145: 1000 format(3(e18.12,2x))
1146: close(44)
1148: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1149: sum = 0.0d+0
1150: do i = 1,mx
1151: xloc = xl0 + (dble(i) * dx)
1152: sum = sum + (x(i) - rexact(xloc,time)) ** 2
1153: enddo
1154: sum = sqrt(sum)
1156: write(*,*)
1157: write(*,*) 'l2norm of the density error is',sum
1158: write(*,*)
1159: endif
1162: return
1163: end
1164: subroutine FormDampit(snes,ctx,x,f,g,y,w, fnorm,ynorm,gnorm,flag,ierr)
1165: ! ---------------------------------------------------------------------
1166: !
1167: ! FormDampit - damps the Newton update
1168: !
1169: ! Input Parameters:
1170: ! snes - the SNES context
1171: ! x - current iterate
1172: ! f - residual evaluated at x
1173: ! y - search direction (containes new iterate on output)
1174: ! w - work vector
1175: ! fnorm - 2-norm of f
1176: !
1177: ! In this example the application context is a Fortran integer array:
1178: ! ctx(1) = shell preconditioner pressure matrix contex
1179: ! ctx(2) = semi implicit pressure matrix
1180: ! ctx(4) = xold - old time values need for time advancement
1181: ! ctx(5) = mx - number of control volumes
1182: ! ctx(6) = N - total number of unknowns
1183: !
1184: ! Output Parameter:
1185: ! g - residual evaluated at new iterate y
1186: ! y - new iterate (contains search direction on input
1187: ! gnorm - 2-norm of g
1188: ! ynorm - 2-norm of search length
1189: ! flag - set to 0 if the line search succeeds; -1 on failure
1190: !
1191: ! Notes:
1192: ! This routine serves as a wrapper for the lower-level routine
1193: ! "ApplicationDampit", where the actual computations are
1194: ! done using the standard Fortran style of treating the local
1195: ! vector data as a multidimensional array over the local mesh.
1196: ! This routine merely accesses the local vector data via
1197: ! VecGetArray() and VecRestoreArray().
1198: !
1199: implicit none
1201: #include finclude/petsc.h
1203: ! Input/output variables:
1204: SNES snes
1205: Vec x, f, g, y, w
1206: PetscFortranAddr ctx(*)
1207: PetscScalar fnorm, ynorm, gnorm
1208: integer ierr, flag
1210: ! Common blocks:
1212: #include "ex74fcomd.h"
1214: ! Local variables:
1216: ! Declarations for use with local arrays:
1217: PetscScalar lx_v(0:1), ly_v(0:1), lw_v(0:1)
1218: PetscOffset lx_i, ly_i , lw_i
1220: !
1221: ! set ynorm
1222: !
1223: call VecNorm(y,NORM_2,ynorm,ierr)
1224: !
1225: ! copy x to w
1226: !
1227: call VecCopy(x,w,ierr)
1228: !
1229: ! get pointers to x, y, w
1230: !
1232: call VecGetArray(x,lx_v,lx_i,ierr)
1233: call VecGetArray(y,ly_v,ly_i,ierr)
1234: call VecGetArray(w,lw_v,lw_i,ierr)
1235: !
1236: ! Compute Damping
1237: !
1238: call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1239: !
1240: ! Restore vectors x, y, w
1241: !
1242: call VecRestoreArray(x,lx_v,lx_i,ierr)
1243: call VecRestoreArray(y,ly_v,ly_i,ierr)
1244: call VecRestoreArray(w,lw_v,lw_i,ierr)
1245: !
1246: ! copy w to y
1247: !
1248: call VecCopy(w,y,ierr)
1249: !
1250: ! compute new residual
1251: !
1252: call SNESComputeFunction(snes,y,g,ierr)
1253: call VecNorm(g,NORM_2,gnorm,ierr)
1254: flag = 0
1256: if (debug) then
1257: write(*,*) 'form damp ynorm = ',ynorm
1258: write(*,*) 'gnorm = ',gnorm
1259: endif
1261: return
1262: end
1263: subroutine FormDt(snes,x,ctx,ierr)
1264: ! ---------------------------------------------------------------------
1265: !
1266: ! FormDt - Compute CFL numbers
1267: !
1268: ! Input Parameters:
1269: ! snes - the SNES context
1270: ! x - input vector
1271: !
1272: ! In this example the application context is a Fortran integer array:
1273: ! ctx(1) = shell preconditioner pressure matrix contex
1274: ! ctx(2) = semi implicit pressure matrix
1275: ! ctx(4) = xold - old time values need for time advancement
1276: ! ctx(5) = mx - number of control volumes
1277: ! ctx(6) = N - total number of unknowns
1278: !
1279: !
1280: ! Notes:
1281: ! This routine serves as a wrapper for the lower-level routine
1282: ! "ApplicationDt", where the actual computations are
1283: ! done using the standard Fortran style of treating the local
1284: ! vector data as a multidimensional array over the local mesh.
1285: ! This routine merely accesses the local vector data via
1286: ! VecGetArray() and VecRestoreArray().
1287: !
1288: implicit none
1290: #include finclude/petsc.h
1292: ! Input/output variables:
1293: SNES snes
1294: Vec x
1295: PetscFortranAddr ctx(*)
1296: integer ierr
1298: ! Common blocks:
1300: #include "ex74fcomd.h"
1302: ! Local variables:
1304: ! Declarations for use with local arrays:
1305: PetscScalar lx_v(0:1)
1306: PetscOffset lx_i
1307: PetscScalar lxold_v(0:1)
1308: PetscOffset lxold_i
1310: !
1311: ! get pointers to x, xold
1312: !
1314: call VecGetArray(x,lx_v,lx_i,ierr)
1315: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1316: !
1317: ! Compute function
1318: !
1319: call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1320: !
1321: ! Restore vectors x, xold
1322: !
1323: call VecRestoreArray(x,lx_v,lx_i,ierr)
1324: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1326: return
1327: end
1328: subroutine FormExact(x,ierr)
1329: ! ---------------------------------------------------------------------
1330: !
1331: ! FormExact - Forms exact solution
1332: !
1333: ! Input Parameter:
1334: ! x - vector
1335: !
1336: ! Output Parameters:
1337: ! x - vector
1338: ! ierr - error code
1339: !
1340: ! Notes:
1341: ! This routine serves as a wrapper for the lower-level routine
1342: ! "ApplicationExact", where the actual computations are
1343: ! done using the standard Fortran style of treating the local
1344: ! vector data as a multidimensional array over the local mesh.
1345: ! This routine merely accesses the local vector data via
1346: ! VecGetArray() and VecRestoreArray().
1347: !
1348: implicit none
1350: #include finclude/petsc.h
1352: ! Input/output variables:
1353: Vec x
1354: integer ierr
1356: ! Declarations for use with local arrays:
1357: PetscScalar lx_v(0:1)
1358: PetscOffset lx_i
1360: 0
1362: !
1363: ! get a pointer to x
1364: !
1365: call VecGetArray(x,lx_v,lx_i,ierr)
1366: !
1367: ! Compute initial guess
1368: !
1369: call ApplicationExact(lx_v(lx_i),ierr)
1370: !
1371: ! Restore vector x
1372: !
1373: call VecRestoreArray(x,lx_v,lx_i,ierr)
1375: return
1376: end
1377: subroutine FormFunction(snes,x,f,ctx,ierr)
1378: ! ---------------------------------------------------------------------
1379: !
1380: ! FormFunction - Evaluates nonlinear function, f(x).
1381: !
1382: ! Input Parameters:
1383: ! snes - the SNES context
1384: ! x - input vector
1385: !
1386: ! In this example the application context is a Fortran integer array:
1387: ! ctx(1) = shell preconditioner pressure matrix contex
1388: ! ctx(2) = semi implicit pressure matrix
1389: ! ctx(4) = xold - old time values need for time advancement
1390: ! ctx(5) = mx - number of control volumes
1391: ! ctx(6) = N - total number of unknowns
1392: !
1393: ! Output Parameter:
1394: ! f - vector with newly computed function
1395: !
1396: ! Notes:
1397: ! This routine serves as a wrapper for the lower-level routine
1398: ! "ApplicationFunction", where the actual computations are
1399: ! done using the standard Fortran style of treating the local
1400: ! vector data as a multidimensional array over the local mesh.
1401: ! This routine merely accesses the local vector data via
1402: ! VecGetArray() and VecRestoreArray().
1403: !
1404: implicit none
1406: #include finclude/petsc.h
1408: ! Input/output variables:
1409: SNES snes
1410: Vec x, f
1411: PetscFortranAddr ctx(*)
1412: integer ierr
1414: ! Common blocks:
1416: #include "ex74fcomd.h"
1418: ! Local variables:
1420: ! Declarations for use with local arrays:
1421: PetscScalar lx_v(0:1), lf_v(0:1)
1422: PetscOffset lx_i, lf_i
1423: PetscScalar lxold_v(0:1)
1424: PetscOffset lxold_i
1426: !
1427: ! get pointers to x, f, xold
1428: !
1430: call VecGetArray(x,lx_v,lx_i,ierr)
1431: call VecGetArray(f,lf_v,lf_i,ierr)
1432: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1433: !
1434: ! Compute function
1435: !
1436: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i), lxold_v(lxold_i),ierr)
1437: !
1438: ! Restore vectors x, f, xold
1439: !
1440: call VecRestoreArray(x,lx_v,lx_i,ierr)
1441: call VecRestoreArray(f,lf_v,lf_i,ierr)
1442: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1443: !
1444: ! something to do with profiling
1445: !
1446: call PetscLogFlops(11*mx,ierr)
1448: return
1449: end
1450: subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1451: ! ---------------------------------------------------------------------
1452: !
1453: ! FormGraph - Forms Graphics output
1454: !
1455: ! Input Parameter:
1456: ! x - vector
1457: ! time - scalar
1458: !
1459: ! Output Parameters:
1460: ! ierr - error code
1461: !
1462: ! Notes:
1463: ! This routine serves as a wrapper for the lower-level routine
1464: ! "ApplicationXmgr", where the actual computations are
1465: ! done using the standard Fortran style of treating the local
1466: ! vector data as a multidimensional array over the local mesh.
1467: ! This routine merely accesses the local vector data via
1468: ! VecGetArray() and VecRestoreArray().
1469: !
1470: implicit none
1472: #include finclude/petsc.h
1474: #include "ex74fcomd.h"
1475: #include "ex74ftube.h"
1477: ! Input/output variables:
1478: Vec x
1479: integer ierr
1481: ! Declarations for use with local arrays:
1482: IS rfrom, rto, rufrom, ruto, efrom, eto
1483: Vec rval
1484: Vec uval
1485: Vec ruval
1486: Vec eval
1487: Vec seval
1488: Vec pval
1489: Vec kval
1490: Vec tval
1491: Vec steval
1492: VecScatter scatter
1493: PetscViewer view0, view1, view2, view3, view4
1494: double precision minus1, l2err, gm1, csubvi
1497: csubvi = 1.0d+0 / csubv
1498: gm1 = gamma - 1.0d+0
1499: 0
1500: minus1 = -1.0d+0
1501: !
1502: ! graphics vectors
1503: !
1504: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,rval,ierr)
1505: call VecSetFromOptions(rval,ierr)
1506: call VecDuplicate(rval,uval,ierr)
1507: call VecDuplicate(rval,ruval,ierr)
1508: call VecDuplicate(rval,eval,ierr)
1509: call VecDuplicate(rval,seval,ierr)
1510: call VecDuplicate(rval,pval,ierr)
1511: call VecDuplicate(rval,kval,ierr)
1512: call VecDuplicate(rval,tval,ierr)
1513: call VecDuplicate(rval,steval,ierr)
1514: !
1515: ! create index sets for vector scatters
1516: !
1517: call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1518: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, rto, ierr)
1519: call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1520: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, ruto, ierr)
1521: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1522: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, eto, ierr)
1524: !
1525: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1526: !
1527: !
1528: ! load rval from x
1529: !
1530: call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1531: call VecScatterBegin(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1532: call VecScatterEnd(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1533: call VecScatterDestroy(scatter,ierr)
1534: !
1535: ! plot rval vector
1536: !
1537: call VecView(rval,view0,ierr)
1538: !
1539: ! make xmgr plot of rval
1540: !
1541: call FormXmgr(rval,1,ierr)
1542: !
1543: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1544: !
1545: !
1546: ! load eval from x
1547: !
1548: call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1549: call VecScatterBegin(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1550: call VecScatterEnd(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1551: call VecScatterDestroy(scatter,ierr)
1552: !
1553: ! plot eval vector
1554: !
1555: call VecView(eval,view2,ierr)
1556: !
1557: ! make xmgr plot of eval
1558: !
1559: call FormXmgr(eval,3,ierr)
1560: !
1561: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1562: !
1564: !
1565: ! load ruval from x
1566: !
1567: call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1568: call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1569: call VecScatterEnd(scatter,x,ruval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1570: call VecScatterDestroy(scatter,ierr)
1571: !
1572: ! create u = ru / r
1573: !
1574: call VecPointwiseDivide(uval,ruval,rval,ierr)
1575: !
1576: ! plot uval vector
1577: !
1578: call VecView(uval,view1,ierr)
1579: !
1580: ! make xmgr plot of uval
1581: !
1582: call FormXmgr(uval,2,ierr)
1584: !
1585: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1586: !
1588: call VecPointwiseMult(kval,uval,uval,ierr)
1589: call VecScale(kval,0.5d+0,ierr)
1591: call VecPointwiseDivide(steval,eval,rval,ierr)
1592: call VecWAXPY(seval,-1.0d+0,kval,steval,ierr)
1594: call VecCopy(seval,tval,ierr)
1595: call VecScale(tval,csubvi,ierr)
1597: !
1598: ! plot tval vector
1599: !
1600: call VecView(tval,view3,ierr)
1601: !
1602: ! make xmgr plot of tval
1603: !
1604: call FormXmgr(tval,4,ierr)
1606: !
1607: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1608: !
1610: call VecPointwiseMult(pval,rval,seval,ierr)
1611: call VecScale(pval,gm1,ierr)
1612: !
1613: ! plot pval vector
1614: !
1615: call VecView(pval,view4,ierr)
1616: !
1617: ! make xmgr plot of pval
1618: !
1619: call FormXmgr(pval,5,ierr)
1620: !
1621: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1622: !
1628: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1629: ! Free work space. All PETSc objects should be destroyed when they
1630: ! are no longer needed.
1631: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1633: call VecDestroy(rval, ierr)
1634: call VecDestroy(uval, ierr)
1635: call VecDestroy(ruval,ierr)
1636: call VecDestroy(eval, ierr)
1637: call VecDestroy(seval, ierr)
1638: call VecDestroy(pval, ierr)
1639: call VecDestroy(kval, ierr)
1640: call VecDestroy(tval, ierr)
1641: call VecDestroy(steval, ierr)
1643: call ISDestroy(rfrom, ierr)
1644: call ISDestroy(rto, ierr)
1646: call ISDestroy(rufrom,ierr)
1647: call ISDestroy(ruto, ierr)
1649: call ISDestroy(efrom, ierr)
1650: call ISDestroy(eto, ierr)
1653: return
1654: end
1655: subroutine FormInitialGuess(x,ierr)
1656: ! ---------------------------------------------------------------------
1657: !
1658: ! FormInitialGuess - Forms initial approximation.
1659: !
1660: ! Input Parameter:
1661: ! x - vector
1662: !
1663: ! Output Parameters:
1664: ! x - vector
1665: ! ierr - error code
1666: !
1667: ! Notes:
1668: ! This routine serves as a wrapper for the lower-level routine
1669: ! "ApplicationInitialGuess", where the actual computations are
1670: ! done using the standard Fortran style of treating the local
1671: ! vector data as a multidimensional array over the local mesh.
1672: ! This routine merely accesses the local vector data via
1673: ! VecGetArray() and VecRestoreArray().
1674: !
1675: implicit none
1677: #include finclude/petsc.h
1679: ! Input/output variables:
1680: Vec x
1681: integer ierr
1683: ! Declarations for use with local arrays:
1684: PetscScalar lx_v(0:1)
1685: PetscOffset lx_i
1687: 0
1689: !
1690: ! get a pointer to x
1691: !
1692: call VecGetArray(x,lx_v,lx_i,ierr)
1693: !
1694: ! Compute initial guess
1695: !
1696: call ApplicationInitialGuess(lx_v(lx_i),ierr)
1697: !
1698: ! Restore vector x
1699: !
1700: call VecRestoreArray(x,lx_v,lx_i,ierr)
1702: return
1703: end
1704: subroutine FormXmgr(x,ivar,ierr)
1705: ! ---------------------------------------------------------------------
1706: !
1707: ! FormXmgr - Forms Xmgr output
1708: !
1709: ! Input Parameter:
1710: ! x - vector
1711: !
1712: ! Output Parameters:
1713: ! x - vector
1714: ! ierr - error code
1715: !
1716: ! Notes:
1717: ! This routine serves as a wrapper for the lower-level routine
1718: ! "ApplicationXmgr", where the actual computations are
1719: ! done using the standard Fortran style of treating the local
1720: ! vector data as a multidimensional array over the local mesh.
1721: ! This routine merely accesses the local vector data via
1722: ! VecGetArray() and VecRestoreArray().
1723: !
1724: implicit none
1726: #include finclude/petsc.h
1728: ! Input/output variables:
1729: Vec x
1730: integer ivar,ierr
1732: ! Declarations for use with local arrays:
1733: PetscScalar lx_v(0:1)
1734: PetscOffset lx_i
1736: 0
1738: !
1739: ! get a pointer to x
1740: !
1741: call VecGetArray(x,lx_v,lx_i,ierr)
1742: !
1743: ! make the graph
1744: !
1745: call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1746: !
1747: ! Restore vector x
1748: !
1749: call VecRestoreArray(x,lx_v,lx_i,ierr)
1751: return
1752: end
1753: subroutine PCRadApply(pc,ctx,x,y,ierr)
1754: ! -------------------------------------------------------------------
1755: !
1756: ! PCRadApply - This routine demonstrates the use of a
1757: ! user-provided preconditioner.
1758: !
1759: ! Input Parameters:
1760: ! dummy - optional user-defined context, not used here
1761: ! x - input vector
1762: ! In this example the shell preconditioner application context
1763: ! is a Fortran integer array:
1764: ! ctx(1) = shell preconditioner pressure matrix contex
1765: ! ctx(2) = semi implicit pressure matrix
1766: ! ctx(4) = xold - old time values need for time advancement
1767: ! ctx(5) = mx - number of control volumes
1768: ! ctx(6) = N - total number of unknowns
1769: !
1770: ! Output Parameters:
1771: ! y - preconditioned vector
1772: ! ierr - error code (nonzero if error has been detected)
1773: !
1774: ! Notes:
1775: ! This code implements the Jacobi preconditioner plus the
1776: ! SOR preconditioner
1777: !
1779: implicit none
1781: #include finclude/petsc.h
1783: #include "ex74fcomd.h"
1784: !
1785: ! Input
1786: !
1787: PC pc
1788: PetscFortranAddr ctx(*)
1789: Vec x, y
1790: integer ierr
1791: !
1792: ! Local
1793: !
1794: IS defrom, deto
1795: Vec de, rese
1796: VecScatter scatter
1797: PetscScalar lde_v(0:1),lrese_v(0:1)
1798: PetscOffset lde_i, lrese_i
1799: !
1800: ! Identity preconditioner
1801: !
1802: call VecCopy(x,y,ierr)
1803: !
1804: ! if kappa0 not equal to zero then precondition the radiation diffusion
1805: !
1806: if (kappa0 .ne. 0.0d+0) then
1807:
1809: !
1810: ! Create needed vectors
1811: !
1812: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,de,ierr)
1813: call VecSetFromOptions(de,ierr)
1814: call VecDuplicate(de,rese,ierr)
1815: !
1816: ! create index sets for scatters
1817: !
1818: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1819: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1820: !
1821: ! load rese from x
1822: !
1823: call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1824: call VecScatterBegin(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1825: call VecScatterEnd(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1826: call VecScatterDestroy(scatter,ierr)
1827: !
1828: ! apply preconditioner
1829: !
1830: call PCApply(ctx(1),rese,de,ierr)
1832: if (debug) then
1833: write(*,*) 'PCRadApply dh is'
1834: call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1835: endif
1836: !
1837: ! load de into y
1838: !
1839: call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1840: call VecScatterBegin(scatter,de,y,INSERT_VALUES, SCATTER_FORWARD,ierr)
1841: call VecScatterEnd(scatter,de,y,INSERT_VALUES,SCATTER_FORWARD,ierr)
1842: call VecScatterDestroy(scatter,ierr)
1844: if (debug) then
1845: write(*,*) 'PCRadApply y is'
1846: call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1847: endif
1849: call VecDestroy(de,ierr)
1850: call VecDestroy(rese,ierr)
1852: call ISDestroy(defrom,ierr)
1853: call ISDestroy(deto,ierr)
1855: endif
1858: return
1859: end
1860: subroutine PCRadSetUp(pc,ctx,ierr)
1861: !
1862: ! PCRadSetUp - This routine sets up a user-defined
1863: ! preconditioner context.
1864: !
1865: ! Input Parameters:
1866: ! In this example the shell preconditioner application context
1867: ! is a Fortran integer array:
1868: ! ctx(1) = shell preconditioner pressure matrix contex
1869: ! ctx(2) = semi implicit pressure matrix
1870: ! ctx(4) = xold - old time values need for time advancement
1871: ! ctx(5) = mx - number of control volumes
1872: ! ctx(6) = N - total number of unknowns
1873: !
1874: ! Output Parameter:
1875: ! ierr - error code (nonzero if error has been detected)
1876: !
1877: ! Notes:
1878: ! In this example, we define the shell preconditioner to be Jacobi
1879: ! method. Thus, here we create a work vector for storing the reciprocal
1880: ! of the diagonal of the preconditioner matrix; this vector is then
1881: ! used within the routine PCRadApply().
1882: !
1884: implicit none
1886: #include finclude/petsc.h
1888: #include "ex74fcomd.h"
1889: !
1890: ! Input
1891: !
1892: PC pc
1893: PetscFortranAddr ctx(*)
1894: integer ierr
1895: !
1896: ! Local
1897: !
1898: Vec eold
1899:
1900: PetscScalar le_v(0:1)
1901: PetscOffset le_i
1902:
1903: !
1904: ! create vector
1905: !
1906: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,eold,ierr)
1907: call VecSetFromOptions(eold,ierr)
1908: !
1909: ! set up the matrix based on xold
1910: !
1911: call Setmat(ctx,ierr)
1912: !
1913: ! set up the preconditioner
1914: !
1915: call PCDestroy(ctx(1),ierr)
1916: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
1917: !VAM call PCSetType(ctx(1),PCJACOBI,ierr)
1918: call PCSetType(ctx(1),PCLU,ierr)
1919: ! call PCSetVector(ctx(1),eold,ierr)
1920: call PCSetOperators(ctx(1),ctx(2),ctx(2), ierr)
1921: call PCSetUp(ctx(1),ierr)
1923: call VecDestroy(eold,ierr)
1926: return
1927: end
1928: subroutine Setmat(ctx,ierr)
1930: implicit none
1932: #include finclude/petsc.h
1934: ! Common blocks:
1935: #include "ex74fcomd.h"
1936: #include "ex74ftube.h"
1938: ! Input/output variables:
1939: PetscFortranAddr ctx(*)
1940: integer ierr
1942: ! Local variables:
1943: PetscScalar lx_v(0:1)
1944: PetscOffset lx_i
1946: double precision xmult, himh, hiph, diag, upper, lower
1947: double precision hi, hip1, him1
1948: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw
1949: double precision see, sep, sew, seef, sewf, tef, twf, ref, rwf, kef, kwf, xmulte, xmultw
1950: !
1951: integer im, nx, ny
1952: !
1953: ! get pointers to xold
1954: !
1955: call VecGetArray(ctx(4),lx_v,lx_i,ierr)
1956:
1958: !
1959: !############################
1960: !
1961: ! loop over all cells begin
1962: !
1963: !############################
1964: !
1965: do im = 1,mx
1966: !
1967: ! set scalars
1968: !
1969: call Setpbc(im,lx_v(lx_i), rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw)
1970: !
1971: ! set diffusion coefficients
1972: !
1973: see = (erge/rhoe) - (0.5d+0 * ue * ue)
1974: sep = (ergp/rhop) - (0.5d+0 * up * up)
1975: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
1977: seef = 0.5d+0 * (see + sep)
1978: sewf = 0.5d+0 * (sew + sep)
1980: tef = seef / csubv
1981: twf = sewf / csubv
1983: ref = 0.5d+0 * (rhoe + rhop)
1984: rwf = 0.5d+0 * (rhow + rhop)
1986: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
1987: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
1989: if (wilson) then
1990: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
1991: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
1992: endif
1993: !
1994: ! set coefficients
1995: !
1996: xmult = dt / (dx * dx * csubv)
1998: xmulte = xmult * kef
1999: xmultw = xmult * kwf
2001: upper = -(xmulte / rhoe)
2002: lower = -(xmultw / rhow)
2004: diag = 1.0d+0 + ( (xmulte + xmultw) / rhop )
2006: !
2007: ! load coefficients into the matrix
2008: !
2009: call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)
2011: if (im .eq. 1) then
2012: call MatSetValues(ctx(2),1,im-1,1,im ,upper, INSERT_VALUES,ierr)
2013: elseif (im .eq. mx) then
2014: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,INSERT_VALUES,ierr)
2015: else
2016: call MatSetValues(ctx(2),1,im-1,1,im ,upper,INSERT_VALUES,ierr)
2017: call MatSetValues(ctx(2),1,im-1,1,im-2,lower, INSERT_VALUES,ierr)
2018: endif
2021: enddo
2022: !
2023: !############################
2024: !
2025: ! loop over all cells end
2026: !
2027: !############################
2028: !
2029:
2030: !
2031: ! final load of matrix
2032: !
2033: call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2034: call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2036: if (debug) then
2037: call MatGetSize(ctx(2),nx,ny,ierr)
2038: write(*,*) 'in setup nx = ',nx,' ny = ',ny
2039: call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2040: endif
2042: call VecRestoreArray (ctx(4),lx_v,lx_i,ierr)
2046: return
2047: end
2048: subroutine Setpbc(i,x, rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw)
2050: implicit none
2052: ! Common blocks:
2053: #include "ex74fcomd.h"
2055: ! Input/output variables:
2056: PetscScalar x(mx*neq)
2057: integer i
2058: double precision rhoee, rhoe, rhop, rhow, rhoww
2059: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2060: double precision ergee, erge, ergp, ergw, ergww
2061: double precision vele, velp, velw
2063: ! Local variables:
2064: integer jr, jru, je
2066: !
2067: ! set pointers
2068: !
2069: jr = (neq*i) - 2
2070: jru = (neq*i) - 1
2071: je = (neq*i)
2073: if (debug) then
2074: write(*,*)
2075: write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2076: write(*,*)
2077: endif
2078:
2079: if (i .eq. 1) then
2081: rhoee = x(jr+(2*neq))
2082: rhoe = x(jr+neq)
2083: rhop = x(jr)
2084: rhow = x(jr)
2085: rhoww = x(jr)
2087: rhouee = x(jru+(2*neq))
2088: rhoue = x(jru+neq)
2089: rhoup = x(jru)
2090: rhouw = x(jru)
2091: rhouww = x(jru)
2093: ergee = x(je+(2*neq))
2094: erge = x(je+neq)
2095: ergp = x(je)
2096: ergw = x(je)
2097: ergww = x(je)
2099: velw = 0.0d+0
2100: velp = rhoup/rhop
2101: vele = rhoue/rhoe
2103: elseif (i .eq. 2) then
2105: rhoee = x(jr+(2*neq))
2106: rhoe = x(jr+neq)
2107: rhop = x(jr)
2108: rhow = x(jr-neq)
2109: rhoww = x(jr-neq)
2111: rhouee = x(jru+(2*neq))
2112: rhoue = x(jru+neq)
2113: rhoup = x(jru)
2114: rhouw = x(jru-neq)
2115: rhouww = x(jru-neq)
2117: ergee = x(je+(2*neq))
2118: erge = x(je+neq)
2119: ergp = x(je)
2120: ergw = x(je-neq)
2121: ergww = x(je-neq)
2123: velw = rhouw/rhow
2124: velp = rhoup/rhop
2125: vele = rhoue/rhoe
2127: elseif (i .eq. mx-1) then
2129: rhoee = x(jr+neq)
2130: rhoe = x(jr+neq)
2131: rhop = x(jr)
2132: rhow = x(jr-neq)
2133: rhoww = x(jr-(2*neq))
2135: rhouee = x(jru+neq)
2136: rhoue = x(jru+neq)
2137: rhoup = x(jru)
2138: rhouw = x(jru-neq)
2139: rhouww = x(jru-(2*neq))
2141: ergee = x(je+neq)
2142: erge = x(je+neq)
2143: ergp = x(je)
2144: ergw = x(je-neq)
2145: ergww = x(je-(2*neq))
2147: velw = rhouw/rhow
2148: velp = rhoup/rhop
2149: vele = rhoue/rhoe
2151: elseif (i .eq. mx) then
2153: rhoee = x(jr)
2154: rhoe = x(jr)
2155: rhop = x(jr)
2156: rhow = x(jr-neq)
2157: rhoww = x(jr-(2*neq))
2159: rhouee = x(jru)
2160: rhoue = x(jru)
2161: rhoup = x(jru)
2162: rhouw = x(jru-neq)
2163: rhouww = x(jru-(2*neq))
2165: ergee = x(je)
2166: erge = x(je)
2167: ergp = x(je)
2168: ergw = x(je-neq)
2169: ergww = x(je-(2*neq))
2171: velw = rhouw/rhow
2172: velp = rhoup/rhop
2173: vele = 0.0d+0
2175: else
2177: rhoee = x(jr+(2*neq))
2178: rhoe = x(jr+neq)
2179: rhop = x(jr)
2180: rhow = x(jr-neq)
2181: rhoww = x(jr-(2*neq))
2183: rhouee = x(jru+(2*neq))
2184: rhoue = x(jru+neq)
2185: rhoup = x(jru)
2186: rhouw = x(jru-neq)
2187: rhouww = x(jru-(2*neq))
2189: ergee = x(je+(2*neq))
2190: erge = x(je+neq)
2191: ergp = x(je)
2192: ergw = x(je-neq)
2193: ergww = x(je-(2*neq))
2195: velw = rhouw/rhow
2196: velp = rhoup/rhop
2197: vele = rhoue/rhoe
2199: endif
2201: if (debug) then
2202: write(*,*)
2203: write(*,*) 'in Setpbc ',i,jr,jru,je
2204: write(*,*) 'mx = ',mx
2205: write(*,*)
2206: endif
2209: return
2210: end
2211: subroutine Setpbcn(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw, jbc)
2213: implicit none
2215: ! Common blocks:
2216: #include "ex74fcomd.h"
2218: ! Input/output variables:
2219: integer jbc
2220: double precision rhoee, rhoe, rhop, rhow, rhoww
2221: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2222: double precision ergee, erge, ergp, ergw, ergww
2223: double precision ue, up, uw
2225: ! Local variables:
2227: if (jbc .eq. 1) then
2228: rhoww = rhop
2229: rhow = rhop
2230: rhouww = rhoup
2231: rhouw = rhoup
2232: ergww = ergp
2233: ergw = ergp
2234: uw = 0.0d+0
2235: elseif (jbc .eq. 2) then
2236: rhoww = rhow
2237: rhouww = rhouw
2238: ergww = ergw
2239: uw = rhouw / rhow
2240: else
2241: uw = rhouw / rhow
2242: endif
2244: if (jbc .eq. mx) then
2245: rhoee = rhop
2246: rhoe = rhop
2247: rhouee = rhoup
2248: rhoue = rhoup
2249: ergee = ergp
2250: erge = ergp
2251: ue = 0.0d+0
2252: elseif (jbc .eq. mx-1) then
2253: rhoee = rhoe
2254: rhouee = rhoue
2255: ergee = erge
2256: ue = rhoue / rhoe
2257: else
2258: ue = rhoue / rhoe
2259: endif
2261: up = rhoup / rhop
2263: if (debug) then
2264: write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2265: endif
2268: return
2269: end
2270: double precision function cont(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww,jcont,xold)
2271: !
2272: ! This function computes the residual
2273: ! for the 1-D continuity equation
2274: !
2275: !
2276: implicit none
2278: include 'ex74fcomd.h'
2279: include 'ex74ftube.h'
2280: !
2281: ! input variables
2282: !
2283: double precision rhoee, rhoe, rhop, rhow, rhoww
2284: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2285: double precision ergee, erge, ergp, ergw, ergww
2286: double precision xold(mx*neq)
2287: !
2288: integer jcont
2289: !
2290: ! local variables
2291: !
2292: double precision theta1
2293: integer jr
2294: !
2295: ! new
2296: !
2297: double precision velfw, velfe
2298: double precision vele,velp,velw
2299: double precision fluxe, fluxw
2300: double precision urhoe, urhow
2301: double precision source
2302: !
2303: ! old
2304: !
2305: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2306: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2307: double precision ergoee, ergoe, ergop, ergow, ergoww
2308: double precision teoee, teoe, teop, teow, teoww, uoe, uoee, uop, uow, uoww
2309: double precision velfow, velfoe
2310: double precision veloe,velop,velow
2311: double precision fluxoe, fluxow
2312: double precision urhooe, urhoow
2313: double precision sourceo
2314: !
2315: ! functions
2316: !
2317: double precision godunov2
2318: double precision upwind, fluxlim
2319: !
2320: !
2321: ! ******************************************************************
2322: !
2323: !
2324: !
2325: if (debug) then
2326: write(*,*)
2327: write(*,*) 'in cont',jcont,' ihod = ',ihod
2328: write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2329: write(*,*) 'rhop = ',rhop
2330: write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2331: write(*,*)
2332: endif
2334: jr = (neq*jcont) - 2
2336: !########################
2337: !
2338: ! NEW
2339: !
2340: !########################
2342: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw, jcont)
2344: velfe = 0.5d+0 * (vele + velp)
2345: velfw = 0.5d+0 * (velw + velp)
2347: if (ihod .eq. 1) then
2349: urhoe = upwind(rhop,rhoe,velfe)
2350: urhow = upwind(rhow,rhop,velfw)
2352: elseif (ihod .eq. 2) then
2354: urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2355: urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)
2357: endif
2359: if (ihod .eq. 3) then
2360: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,1)
2361: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,1)
2362: else
2363: fluxe = (dt/dx) * urhoe
2364: fluxw = (dt/dx) * urhow
2365: endif
2368: source = 0.0d+0
2370: !########################
2371: !
2372: ! OLD
2373: !
2374: !########################
2376: call Setpbc(jcont,xold,rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
2378: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jcont)
2380: velfoe = 0.5d+0 * (veloe + velop)
2381: velfow = 0.5d+0 * (velow + velop)
2384: if (ihod .eq. 1) then
2386: urhooe = upwind(rhoop,rhooe,velfoe)
2387: urhoow = upwind(rhoow,rhoop,velfow)
2389: elseif (ihod .eq. 2) then
2391: urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2392: urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)
2394: endif
2396: if (ihod .eq. 3) then
2397: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,1)
2398: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,1)
2399: else
2400: fluxoe = (dt/dx) * urhooe
2401: fluxow = (dt/dx) * urhoow
2402: endif
2404: sourceo = 0.0d+0
2407: !########################
2408: !
2409: ! FUNCTION
2410: !
2411: !########################
2413: theta1 = 1.0d+0 - theta
2414: cont = (rhop - xold(jr)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2415: !VAM
2416: if (probnum .eq. 3) then
2417: cont = 0.0d+0
2418: endif
2419: !VAM
2422: if (debug) then
2423: write(*,*)
2424: write(*,*) 'cont(',jcont,') = ',cont
2425: write(*,*) 'theta = ',theta,'rhop = ',rhop
2426: write(*,*) 'source = ',source,' sourceo = ',sourceo
2427: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2428: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2429: write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2430: write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2431: write(*,*)
2432: endif
2434: return
2435: end
2436: double precision function eexact(x,t)
2438: implicit none
2440: double precision x,t
2441: double precision xot, head, tail, contact, ufan
2442: double precision xpow, grat, urat
2443: double precision uexact
2446: logical debug
2448: include 'ex74ftube.h'
2450: debug = .false.
2453: if (t .le. 0.0d+0) then
2454: if (x .gt. 0.0d+0) then
2455: eexact = e1
2456: else
2457: eexact = e4
2458: endif
2459: else
2461: xot = x/t
2462: head = -a4
2463: tail = v3 - a3
2464: contact = v2
2466: if (xot .lt. head) then
2467: eexact = e4
2468: elseif (xot .gt. sspd) then
2469: eexact = e1
2470: elseif (xot .gt. contact) then
2471: eexact = e2
2472: elseif (xot .gt. tail) then
2473: eexact = e3
2474: else
2475: ufan = uexact(x,t)
2476: grat = (gamma - 1.0d+0) / 2.0d+0
2477: xpow = 2.0d+0
2478: urat = ufan / a4
2479: eexact = e4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
2480: endif
2482: endif
2485: if (debug) then
2486: write(*,*)
2487: write(*,*) 'eexact(',x,',',t,') = ',eexact
2488: write(*,*)
2489: endif
2491: return
2492: end
2493: subroutine eigen(ht,uht)
2494: !23456789012345678901234567890123456789012345678901234567890123456789012
2495: !
2496: ! subroutine eigen
2497: !
2498: ! This subroutine computes the eigen values and eigen vectors
2499: !
2500: !23456789012345678901234567890123456789012345678901234567890123456789012
2503: !#######################################################################
2505: implicit none
2507: include 'ex74fcomd.h'
2509: double precision ht, uht
2511: double precision ut, at, lam1, lam2
2514: !#######################################################################
2516: ut = uht / ht
2517: at = sqrt( ht)
2519: lam1 = ut - at
2520: lam2 = ut + at
2522: eigval(1) = lam1
2523: eigval(2) = lam2
2525: eigvec(1,1) = 1.0d+0
2526: eigvec(2,1) = lam1
2527: eigvec(1,2) = 1.0d+0
2528: eigvec(2,2) = lam2
2530: rinv(1,1) = lam2 / (2.0d+0 * at)
2531: rinv(2,1) = -lam1 / (2.0d+0 * at)
2532: rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2533: rinv(2,2) = 1.0d+0 / (2.0d+0 * at)
2536: return
2537: end
2538: subroutine eigene(r,ru,e,l1,l2,l3)
2539: !23456789012345678901234567890123456789012345678901234567890123456789012
2540: !
2541: ! subroutine eigene
2542: !
2543: ! This subroutine computes the eigen values for the entropy fix
2544: !
2545: !23456789012345678901234567890123456789012345678901234567890123456789012
2548: !#######################################################################
2550: implicit none
2552: include 'ex74ftube.h'
2554: double precision r,ru,e,l1,l2,l3
2556: double precision p,u,a
2558: double precision eos
2559: integer ierr
2561: logical debug
2564: !#######################################################################
2566: debug = .false.
2568: if (debug) then
2569: write(*,*)
2570: write(*,*) 'gamma = ',gamma
2571: write(*,*) 'r,ru,e = ',r,ru,e
2572: write(*,*)
2573: endif
2575: p = eos(r,ru,e)
2576: u = ru/r
2577: if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2578: write(*,*)
2579: write(*,*) 'gamma = ',gamma
2580: write(*,*) 'r = ',r
2581: write(*,*) 'p = ',p
2582: write(*,*)
2583: call PetscFinalize(ierr)
2584: stop
2585: endif
2586: a = sqrt((gamma * p)/r)
2588: if (debug) then
2589: write(*,*)
2590: write(*,*) 'p,u,a = ',p,u,a
2591: write(*,*)
2592: endif
2594: l1 = u - a
2595: l2 = u
2596: l3 = u + a
2598: return
2599: end
2600: double precision function energy(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,jerg,xold)
2601: !
2602: ! This function computes the residual
2603: ! for the 1-D energy equation
2604: !
2605: !
2606: implicit none
2608: include 'ex74fcomd.h'
2609: include 'ex74ftube.h'
2610: !
2611: ! input variables
2612: !
2613: double precision rhoee, rhoe, rhop, rhow, rhoww
2614: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2615: double precision ergee, erge, ergp, ergw, ergww
2616: double precision xold(mx*neq)
2617: !
2618: integer jerg
2619: !
2620: ! local variables
2621: !
2622: double precision theta1
2623: integer je
2624: !
2625: ! new
2626: !
2627: double precision velfw, velfe
2628: double precision vele,velp,velw
2629: double precision fluxe, fluxw
2630: double precision uepe, uepw
2631: double precision ue, up, uw
2632: double precision see, sep, sew
2633: double precision seef, sewf
2634: double precision upe, upw
2635: double precision presse, pressw
2636: double precision source
2637: double precision te, tp, tw
2638: double precision tef, twf, ref, rwf
2639: double precision kef, kwf
2640: double precision hflxe, hflxw
2641: !
2642: ! old
2643: !
2644: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2645: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2646: double precision ergoee, ergoe, ergop, ergow, ergoww
2647: double precision velfow, velfoe
2648: double precision veloe,velop,velow
2649: double precision fluxoe, fluxow
2650: double precision uepoe, uepow
2651: double precision uoe, uop, uow
2652: double precision seoe, seop, seow
2653: double precision seoef, seowf
2654: double precision upoe, upow
2655: double precision pressoe, pressow
2656: double precision sourceo
2657: double precision toe, top, tow
2658: double precision toef, towf, roef, rowf
2659: double precision koef, kowf
2660: double precision hflxoe, hflxow
2661: !
2662: ! functions
2663: !
2664: double precision godunov2, eos
2665: double precision upwind, fluxlim
2667: !
2668: !
2669: ! ******************************************************************
2670: !
2671: !
2672: !
2673: je = (neq*jerg)
2675: !########################
2676: !
2677: ! NEW
2678: !
2679: !########################
2681: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,vele, velp, velw, jerg)
2683: pressw = eos(rhow, rhouw, ergw)
2684: presse = eos(rhoe, rhoue, erge)
2686: uw = rhouw / rhow
2687: up = rhoup / rhop
2688: ue = rhoue / rhoe
2690: upw = uw * pressw
2691: upe = ue * presse
2693: velfe = 0.5d+0 * (vele + velp)
2694: velfw = 0.5d+0 * (velw + velp)
2696: if (ihod .eq. 1) then
2698: uepe = upwind(ergp,erge,velfe)
2699: uepw = upwind(ergw,ergp,velfw)
2701: elseif (ihod .eq. 2) then
2703: uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2704: uepw = fluxlim(ergww,ergw,ergp,erge,velfw)
2706: endif
2708: if (ihod .eq. 3) then
2709: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,3)
2710: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,3)
2711: else
2712: fluxe = (dt/dx) * ( uepe + (0.5d+0*upe) )
2713: fluxw = (dt/dx) * ( uepw + (0.5d+0*upw) )
2714: endif
2715: !
2716: ! radiation
2717: !
2718: if (kappa0 .eq. 0.0d+0) then
2719: source = 0.0d+0
2720: else
2722: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2723: sep = (ergp/rhop) - (0.5d+0 * up * up)
2724: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2726: seef = 0.5d+0 * (see + sep)
2727: sewf = 0.5d+0 * (sew + sep)
2729: te = see / csubv
2730: tp = sep / csubv
2731: tw = sew / csubv
2733: tef = seef / csubv
2734: twf = sewf / csubv
2736: ref = 0.5d+0 * (rhoe + rhop)
2737: rwf = 0.5d+0 * (rhow + rhop)
2739: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2740: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2742: if (wilson) then
2743: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2744: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2745: endif
2747: if ( debug .and. (kef .gt. 1.0d+10) ) then
2748: write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2749: endif
2750: if ( debug .and. (kwf .gt. 1.0d+10) ) then
2751: write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2752: endif
2754: hflxe = kef * (te - tp) / dx
2755: hflxw = kwf * (tp - tw) / dx
2757: source = (dt/dx) * (hflxe - hflxw)
2759: endif
2761: !########################
2762: !
2763: ! OLD
2764: !
2765: !########################
2767: call Setpbc(jerg,xold, rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
2769: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jerg)
2771: pressow = eos(rhoow, rhouow, ergow)
2772: pressoe = eos(rhooe, rhouoe, ergoe)
2774: uow = rhouow / rhoow
2775: uop = rhouop / rhoop
2776: uoe = rhouoe / rhooe
2778: upow = uow * pressow
2779: upoe = uoe * pressoe
2781: velfoe = 0.5d+0 * (veloe + velop)
2782: velfow = 0.5d+0 * (velow + velop)
2785: if (ihod .eq. 1) then
2787: uepoe = upwind(ergop,ergoe,velfoe)
2788: uepow = upwind(ergow,ergop,velfow)
2790: elseif (ihod .eq. 2) then
2792: uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
2793: uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)
2795: endif
2797: if (ihod .eq. 3) then
2798: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,3)
2799: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,3)
2800: else
2801: fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
2802: fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
2803: endif
2805: !
2806: ! old radiation
2807: !
2808: if (kappa0 .eq. 0.0d+0) then
2809: sourceo = 0.0d+0
2810: else
2812: seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
2813: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
2814: seow = (ergow/rhoow) - (0.5d+0 * uow * uow)
2816: seoef = 0.5d+0 * (seoe + seop)
2817: seowf = 0.5d+0 * (seow + seop)
2819: toe = seoe / csubv
2820: top = seop / csubv
2821: tow = seow / csubv
2823: toef = seoef / csubv
2824: towf = seowf / csubv
2826: roef = 0.5d+0 * (rhooe + rhoop)
2827: rowf = 0.5d+0 * (rhoow + rhoop)
2829: koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
2830: kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)
2832: if (wilson) then
2833: koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
2834: kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
2835: endif
2837: if ( debug .and. (koef .gt. 1.0d+10) ) then
2838: write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
2839: endif
2840: if ( debug .and. (kowf .gt. 1.0d+10) ) then
2841: write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
2842: endif
2844: hflxoe = koef * (toe - top) / dx
2845: hflxow = kowf * (top - tow) / dx
2847: sourceo = (dt/dx) * (hflxoe - hflxow)
2849: endif
2852: !########################
2853: !
2854: ! FUNCTION
2855: !
2856: !########################
2858: !VAM
2859: if (probnum .eq. 3) then
2860: fluxe = 0.0d+0
2861: fluxw = 0.0d+0
2862: fluxoe = 0.0d+0
2863: fluxow = 0.0d+0
2864: endif
2865: !VAM
2867: theta1 = 1.0d+0 - theta
2868: energy = (ergp - xold(je)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2870: if (debug) then
2871: write(*,*)
2872: write(*,*) 'energy(',jerg,') = ',energy
2873: write(*,*)
2874: write(*,*) fluxe,fluxw
2875: write(*,*) fluxoe,fluxow
2876: write(*,*) source,sourceo
2877: write(*,*)
2878: endif
2880: return
2881: end
2882: double precision function eos(r,ru,e)
2884: implicit none
2886: double precision r,ru,e
2888: double precision se, u
2890: integer ierr
2892: logical debug
2894: include "ex74ftube.h"
2896: debug = .false.
2898: if (debug) then
2899: write(*,*)
2900: write(*,*) 'in eos r,ru,e'
2901: write(*,*) r,ru,e
2902: write(*,*)
2903: endif
2905: u = ru/r
2907: se = (e/r) - (0.5d+0 * u * u)
2908: eos = (gamma - 1.0d+0) * r * se
2910: if (eos .lt. 0.0d+0) then
2911: write(*,*)
2912: write(*,*) 'eos = ',eos
2913: write(*,*) 'gamma = ',gamma
2914: write(*,*) 'r = ',r
2915: write(*,*) 'se = ',se
2916: write(*,*) 'e = ',e
2917: write(*,*) 'u = ',u
2918: write(*,*) 'ru = ',ru
2919: call PetscFinalize(ierr)
2920: write(*,*)
2921: stop
2922: endif
2924: if (debug) then
2925: write(*,*)
2926: write(*,*) 'in eos u,se,eos'
2927: write(*,*) u,se,eos
2928: write(*,*)
2929: endif
2932: return
2933: end
2934: subroutine eval2
2936: implicit none
2938: double precision prat, grat, xnum, xdenom
2941: logical debug
2943: include 'ex74ftube.h'
2945: debug = .false.
2947: prat = p2/p1
2948: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
2950: xnum = grat + prat
2951: xdenom = 1.0d+0 + (prat * grat)
2952:
2953: e2 = e1 * prat * (xnum/xdenom)
2954:
2957: if (debug) then
2958: write(*,*)
2959: write(*,*) 'e1 = ',e1
2960: write(*,*) 'e2 = ',e2
2961: endif
2963: return
2964: end
2965: subroutine exact0
2967: implicit none
2969: double precision tol, xn
2970: double precision shockp, fprime
2972: integer maxnewt, niter
2974: logical found, debug
2976: include 'ex74ftube.h'
2978: debug = .false.
2980: tol = 1.0d-10
2982: maxnewt = 40
2983:
2984: a1 = sqrt(gamma*p1/r1)
2985: a4 = sqrt(gamma*p4/r4)
2989: found = .false.
2990: niter = 0
2992: xn = 0.5d+0 * (p1 + p4)
2994: 10 if ( (.not. found) .and. (niter .le. maxnewt) ) then
2996: niter = niter + 1
2998: xn = xn - (shockp(xn) / fprime(xn))
3000: if (debug) then
3001: write(*,*) niter,shockp(xn),xn
3002: endif
3004: if ( abs(shockp(xn)) .lt. tol ) then
3005: found = .true.
3006: endif
3008: goto 10
3010: endif
3012: if (.not. found) then
3014: write(*,*) 'newton failed'
3015: write(*,*) xn,shockp(xn)
3016: stop
3018: endif
3020: p2 = xn
3023: if (debug) then
3024: write(*,*)
3025: write(*,*) 'p1 = ',p1
3026: write(*,*) 'p2 = ',p2
3027: write(*,*) 'p4 = ',p4
3028: write(*,*)
3029: endif
3031: return
3032: end
3033: double precision function flux(r,ru,e,eqn)
3034: !23456789012345678901234567890123456789012345678901234567890123456789012
3035: !
3036: ! function flux
3037: !
3038: ! This function computes the flux at a face
3039: !
3040: !23456789012345678901234567890123456789012345678901234567890123456789012
3043: !#######################################################################
3045: implicit none
3047: include 'ex74fcomd.h'
3048: include 'ex74ftube.h'
3050: double precision r, ru, e
3052: integer eqn
3054: double precision p,u
3056: double precision eos
3059: !#######################################################################
3061: p = eos(r,ru,e)
3062: u = ru/r
3064: if (eqn .eq. 1) then
3065: flux = ru
3066: elseif (eqn .eq. 2) then
3067: flux = (u * ru) + p
3068: else
3069: flux = u * (e + p)
3070: endif
3072: return
3073: end
3074: double precision function fluxlim(fww,fw,fe,fee,vp)
3075: !23456789012345678901234567890123456789012345678901234567890123456789012
3076: !
3077: ! function fluxlim
3078: !
3079: ! this function computes the flux limited quick face value
3080: !
3081: !23456789012345678901234567890123456789012345678901234567890123456789012
3084: !#######################################################################
3086: implicit none
3088: double precision fww, fw, fe, fee, vp
3090: double precision fd, fc, fu
3092: double precision f1, f2, f3, f4, fhod, beta, flc
3094: double precision med, quick
3095:
3096: logical limit
3098: !#######################################################################
3100: limit = .true.
3102: if (vp .gt. 0.0d+0) then
3103: fd = fe
3104: fc = fw
3105: fu = fww
3106: else
3107: fd = fw
3108: fc = fe
3109: fu = fee
3110: endif
3112: fhod = quick(fd,fc,fu)
3114: if (limit) then
3116: beta = 0.25d+0
3117: flc = 4.0d+0
3119: f1 = fc
3120: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3121: f3 = fu + ( flc * (fc - fu) )
3122: f4 = med(f1,f2,f3)
3123: fluxlim = vp * med(f1,f4,fhod)
3125: else
3127: fluxlim = vp * fhod
3129: endif
3131: return
3132: end
3133: double precision function fluxlim2(fww,fw,fe,fee,vp)
3134: !23456789012345678901234567890123456789012345678901234567890123456789012
3135: !
3136: ! function fluxlim2
3137: !
3138: ! this function computes the flux limited quick face value
3139: !
3140: !23456789012345678901234567890123456789012345678901234567890123456789012
3143: !#######################################################################
3145: implicit none
3147: double precision fww, fw, fe, fee, vp
3149: double precision fd, fc, fu
3151: double precision f1, f2, f3, f4, fhod, beta, flc
3153: double precision med, quick
3154:
3155: logical limit, debug
3157: !#######################################################################
3159: debug = .false.
3161: if (debug) then
3162: write(*,*)
3163: write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3164: write(*,*) fee,fe,fw,fww
3165: write(*,*)
3166: endif
3168: limit = .true.
3170: if (vp .gt. 0.0d+0) then
3171: fd = fe
3172: fc = fw
3173: fu = fww
3174: else
3175: fd = fw
3176: fc = fe
3177: fu = fee
3178: endif
3180: fhod = quick(fd,fc,fu)
3182: if (limit) then
3184: beta = 0.25d+0
3185: flc = 4.0d+0
3187: f1 = fc
3188: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3189: f3 = fu + ( flc * (fc - fu) )
3190: f4 = med(f1,f2,f3)
3191: fluxlim2 = med(f1,f4,fhod)
3193: else
3195: fluxlim2 = fhod
3197: endif
3199: return
3200: end
3201: double precision function fprime(x)
3203: implicit none
3205: double precision x, eps
3206: double precision shockp
3208: eps = 1.0d-8
3210: fprime = ( shockp(x+eps) - shockp(x) ) / eps
3212: return
3213: end
3214: double precision function godunov2(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3215: !23456789012345678901234567890123456789012345678901234567890123456789012
3216: !
3217: ! function godunov2
3218: !
3219: ! this function computes the roe/godunov2 face value
3220: !
3221: !23456789012345678901234567890123456789012345678901234567890123456789012
3224: !#######################################################################
3226: implicit none
3228: include 'ex74fcomd.h'
3229: include 'ex74ftube.h'
3231: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3233: integer eqn
3235: double precision rrg, rlg, rurg, rulg, erg, elg
3237: double precision hlle
3240: !#######################################################################
3242: if (gorder .eq. 1) then
3243: rrg = rr
3244: rlg = rl
3245: rurg = rur
3246: rulg = rul
3247: erg = er
3248: elg = el
3249: else
3250: call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,rrg, rlg,rurg, rulg, erg, elg)
3251: endif
3253: !VAM if (ientro .eq. 0) then
3254: !VAM godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3255: !VAM elseif(ientro .eq. 1) then
3256: !VAM godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3257: !VAM else
3258: godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3259: !VAM endif
3262: return
3263: end
3264: double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3265: !23456789012345678901234567890123456789012345678901234567890123456789012
3266: !
3267: ! function hlle
3268: !
3269: ! this function computes the roe/hlle face value
3270: !
3271: !23456789012345678901234567890123456789012345678901234567890123456789012
3274: !#######################################################################
3276: implicit none
3278: include 'ex74fcomd.h'
3279: include 'ex74ftube.h'
3281: double precision rrg,rlg,rurg,rulg,erg,elg
3282: integer eqn
3284: double precision laml1, laml2, laml3
3285: double precision lamr1, lamr2, lamr3
3286: double precision sl, sr
3289: double precision flux
3291: integer i, j, ispeed
3292:
3294: !#######################################################################
3296: ispeed = 1
3298: do i = 1,neq
3299: fr(i) = flux(rrg,rurg,erg,i)
3300: fl(i) = flux(rlg,rulg,elg,i)
3301: enddo
3303: deltau(1) = rrg - rlg
3304: deltau(2) = rurg - rulg
3305: deltau(3) = erg - elg
3307: !VAM call roestat(uhl,uhr,hl,hr,ht,uht)
3309: !VAM call eigene(ht,uht,lamt1, lamt2)
3310: call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3311: call eigene(rlg,rulg,elg,laml1,laml2,laml3)
3313: !VAM if (ispeed .eq. 1) then
3314: !VAM sl = min(laml1,lamt1)
3315: !VAM sr = max(lamt2,lamr2)
3316: !VAM else
3317: sl = min(laml1,lamr1)
3318: sr = max(laml3,lamr3)
3319: !VAM endif
3322: do i = 1,neq
3323: froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )/(sr-sl)
3324: enddo
3326: hlle = froe(eqn)
3329: return
3330: end
3331: double precision function med(x1,x2,x3)
3332: !23456789012345678901234567890123456789012345678901234567890123456789012
3333: !
3334: ! function med
3335: !
3336: ! this function computes the median of three numbers
3337: !
3338: !23456789012345678901234567890123456789012345678901234567890123456789012
3341: !#######################################################################
3343: implicit none
3345: double precision x1, x2, x3
3346: double precision xhi, xlo
3348: !#######################################################################
3350: xhi = max(x1,x2,x3)
3351: xlo = min(x1,x2,x3)
3353: med = x1 + x2 + x3 - xhi - xlo
3355: return
3356: end
3357: double precision function mom(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,jmom,xold)
3358: !
3359: ! This function computes the residual
3360: ! for the 1-D momentum equation
3361: !
3362: !
3363: implicit none
3365: include 'ex74fcomd.h'
3366: include 'ex74ftube.h'
3367: !
3368: ! input variables
3369: !
3370: double precision rhoee, rhoe, rhop, rhow, rhoww
3371: double precision rhouee, rhoue, rhoup, rhouw, rhouww
3372: double precision ergee, erge, ergp, ergw, ergww
3373: double precision xold(mx*neq)
3374: !
3375: integer jmom
3376: !
3377: ! local variables
3378: !
3379: double precision theta1
3380: integer jru
3381: !
3382: ! new
3383: !
3384: double precision velfw, velfe
3385: double precision vele,velp,velw
3386: double precision fluxe, fluxw
3387: double precision uurhoe, uurhow
3388: double precision pressee, presse, pressp,pressw, pressww
3389: double precision rupee, rupe, rupp, rupw, rupww
3390: double precision uee, ue, up, uw, uww
3391: double precision source
3392: !
3393: ! old
3394: !
3395: double precision rhooee, rhooe, rhoop, rhoow, rhooww
3396: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3397: double precision ergoee, ergoe, ergop, ergow, ergoww
3398: double precision velfow, velfoe
3399: double precision veloe,velop,velow
3400: double precision fluxoe, fluxow
3401: double precision uurhooe, uurhoow
3402: double precision pressoee, pressoe, pressop, pressow, pressoww
3403: double precision rupoee, rupoe, rupop, rupow, rupoww
3404: double precision uoee, uoe, uop, uow, uoww
3405: double precision sourceo
3407: double precision eps
3408: !
3409: ! functions
3410: !
3411: double precision godunov2, eos
3412: double precision upwind, fluxlim
3413: !
3414: !
3415: ! ******************************************************************
3416: !
3417: !
3418: eps = 1.0d-32
3419: !
3420: jru = (neq*jmom) - 1
3422: !########################
3423: !
3424: ! NEW
3425: !
3426: !########################
3428: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,vele, velp, velw, jmom)
3430: presse = eos(rhoe, rhoue, erge )
3431: pressw = eos(rhow, rhouw, ergw )
3433: velfe = 0.5d+0 * (vele + velp)
3434: velfw = 0.5d+0 * (velw + velp)
3436: if (ihod .eq. 1) then
3438: uurhoe = upwind(rhoup,rhoue,velfe)
3439: uurhow = upwind(rhouw,rhoup,velfw)
3441: elseif (ihod .eq. 2) then
3443: uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3444: uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)
3446: endif
3448: if (ihod .eq. 3) then
3449: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,2)
3450: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,2)
3451: else
3452: fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3453: fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3454: endif
3457: source = 0.0d+0
3459: !########################
3460: !
3461: ! OLD
3462: !
3463: !########################
3465: call Setpbc(jmom,xold,rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
3467: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jmom)
3469: pressoe = eos(rhooe, rhouoe, ergoe)
3470: pressow = eos(rhoow, rhouow, ergow)
3472: velfoe = 0.5d+0 * (veloe + velop)
3473: velfow = 0.5d+0 * (velow + velop)
3475: if (ihod .eq. 1) then
3477: uurhooe = upwind(rhouop,rhouoe,velfoe)
3478: uurhoow = upwind(rhouow,rhouop,velfow)
3480: elseif (ihod .eq. 2) then
3482: uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3483: uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)
3485: endif
3487: if (ihod .eq. 3) then
3488: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee,ergow, ergop, ergoe, ergoee,2)
3489: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,2)
3490: else
3491: fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3492: fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3493: endif
3495: sourceo = 0.0d+0
3498: !########################
3499: !
3500: ! FUNCTION
3501: !
3502: !########################
3504: theta1 = 1.0d+0 - theta
3505: mom = (rhoup - xold(jru)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3506: !VAM
3507: if (probnum .eq. 3) then
3508: mom = 0.0d+0
3509: endif
3510: !VAM
3511: if (debug) then
3512: write(*,*)
3513: write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3514: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3515: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3516: write(*,*) 'presse = ',presse,'pressw = ',pressw
3517: write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3518: write(*,*)
3519: endif
3521: return
3522: end
3523: double precision function quick(fd, fc, fu)
3524: !23456789012345678901234567890123456789012345678901234567890123456789012
3525: !
3526: ! function quick
3527: !
3528: ! this function computes the quick face value
3529: !
3530: !23456789012345678901234567890123456789012345678901234567890123456789012
3533: !#######################################################################
3535: implicit none
3537: double precision fd, fc, fu
3539: !#######################################################################
3541: quick = ( (3.0d+0 * fd) + (6.0d+0 * fc) - fu ) / 8.0d+0
3543: return
3544: end
3545: double precision function rexact(x,t)
3547: implicit none
3549: double precision x,t
3550: double precision xot, head, tail, contact, ufan
3551: double precision xpow, grat, urat
3552: double precision uexact
3555: logical debug
3557: include 'ex74ftube.h'
3559: debug = .false.
3562: if (t .le. 0.0d+0) then
3563: if (x .gt. 0.0d+0) then
3564: rexact = r1
3565: else
3566: rexact = r4
3567: endif
3568: else
3570: xot = x/t
3571: head = -a4
3572: tail = v3 - a3
3573: contact = v2
3575: if (xot .lt. head) then
3576: rexact = r4
3577: elseif (xot .gt. sspd) then
3578: rexact = r1
3579: elseif (xot .gt. contact) then
3580: rexact = r2
3581: elseif (xot .gt. tail) then
3582: rexact = r3
3583: else
3584: ufan = uexact(x,t)
3585: grat = (gamma - 1.0d+0) / 2.0d+0
3586: xpow = 1.0d+0 / grat
3587: urat = ufan / a4
3588: rexact = r4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
3589: endif
3591: endif
3594: if (debug) then
3595: write(*,*)
3596: write(*,*) 'rexact(',x,',',t,') = ',rexact
3597: write(*,*)
3598: endif
3600: return
3601: end
3602: subroutine roestat(uhl, uhr, hl,hr,ht,uht)
3603: !23456789012345678901234567890123456789012345678901234567890123456789012
3604: !
3605: ! subroutine roestat
3606: !
3607: ! This subroutine computes the roe state at a face
3608: !
3609: !23456789012345678901234567890123456789012345678901234567890123456789012
3612: !#######################################################################
3614: implicit none
3616: include 'ex74fcomd.h'
3618: double precision uhl, uhr, hl, hr, ht, uht
3620: double precision ul, ur, shl, shr, xnum, xdenom
3621:
3624: !#######################################################################
3626: ul = uhl / hl
3627: ur = uhr / hr
3629: shl = sqrt(hl)
3630: shr = sqrt(hr)
3632: xnum = (shl * ul) + (shr * ur)
3633: xdenom = shl + shr
3635: ht = 0.5d+0 * (hl + hr)
3636: uht = ht * ( xnum / xdenom )
3638: return
3639: end
3640: subroutine rval2
3642: implicit none
3644: double precision prat, grat, xnum, xdenom
3647: logical debug
3649: include 'ex74ftube.h'
3651: debug = .false.
3653: prat = p2/p1
3654: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
3656: xnum = 1.0d+0 + (grat * prat)
3657: xdenom = grat + prat
3658:
3659: r2 = r1 * (xnum/xdenom)
3660:
3663: if (debug) then
3664: write(*,*)
3665: write(*,*) 'r1 = ',r1
3666: write(*,*) 'r2 = ',r2
3667: endif
3669: return
3670: end
3671: subroutine secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err, rrg, rlg,rurg, rulg, erg, elg)
3672: !23456789012345678901234567890123456789012345678901234567890123456789012
3673: !
3674: ! subroutine secondq
3675: !
3676: ! this subroutine computes the second order (based on quick) left
3677: ! and right states for the godunov solver.
3678: !
3679: !23456789012345678901234567890123456789012345678901234567890123456789012
3682: !#######################################################################
3684: implicit none
3686: include 'ex74fcomd.h'
3688: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3689: double precision rrg, rlg,rurg, rulg, erg, elg
3693: double precision veld, ull,ul,ur,urr, ulg, urg
3695: double precision fluxlim2
3698: !#######################################################################
3700: !
3701: ! compute the velocities
3702: !
3703: ull = rull/rll
3704: ul = rul /rl
3705: ur = rur /rr
3706: urr = rurr/rrr
3708: !
3709: ! compute the left state first
3710: !
3711: veld = 1.0d+0
3713: rlg = fluxlim2(rll,rl,rr,rrr,veld)
3714: ulg = fluxlim2(ull,ul,ur,urr,veld)
3715: rulg = rlg * ulg
3716: elg = fluxlim2(ell,el,er,err,veld)
3717: !
3718: ! now compute the right state
3719: !
3720: veld = -1.0d+0
3722: rrg = fluxlim2(rll,rl,rr,rrr,veld)
3723: urg = fluxlim2(ull,ul,ur,urr,veld)
3724: rurg = rrg * urg
3725: erg = fluxlim2(ell,el,er,err,veld)
3729: return
3730: end
3731: double precision function shockp(x)
3733: implicit none
3735: double precision x
3736: double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
3737: logical debug
3739: include 'ex74ftube.h'
3741: debug = .false.
3744: if (debug) then
3745: write(*,*)
3746: write(*,*) 'gamma = ',gamma
3747: write(*,*) 'a1 = ',a1
3748: write(*,*) 'a4 = ',a4
3749: write(*,*) 'p1 = ',p1
3750: write(*,*) 'p2 = ',x
3751: write(*,*)
3752: endif
3754: xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
3755: xdenom = sqrt ( 2.0d+0 * gamma * ( (2.0d+0*gamma) + (gamma + 1.0d+0) * ((x/p1) - 1) ) )
3756: xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)
3758: shockp = (x/p1)*((1.0d+0-(xnum/xdenom))**xpow) - (p4/p1)
3761: if (debug) then
3762: write(*,*)
3763: write(*,*) 'xnum = ',xnum
3764: write(*,*) 'gamma = ',gamma
3765: write(*,*) 'a1 = ',a1
3766: write(*,*) 'a4 = ',a4
3767: write(*,*) 'p1 = ',p1
3768: write(*,*) 'xdenom = ',xdenom
3769: write(*,*) 'xpow = ',xpow
3770: write(*,*) 'shockp = ',shockp
3771: write(*,*) 'p2 = ',x
3772: write(*,*)
3773: endif
3775: return
3776: end
3777: double precision function uexact(x,t)
3779: implicit none
3781: double precision x,t
3782: double precision xot, head, tail
3785: logical debug
3787: include 'ex74ftube.h'
3789: debug = .false.
3791: if (debug) then
3792: write(*,*)
3793: write(*,*) 't = ',t
3794: write(*,*) 'x = ',x
3795: write(*,*) 'a4 = ',a4
3796: write(*,*) 'v3 = ',v3
3797: write(*,*) 'a3 = ',a3
3798: write(*,*)
3799: endif
3801: if (t .le. 0.0d+0) then
3802: uexact = 0.0d+0
3803: else
3805: xot = x/t
3806: head = -a4
3807: tail = v3 - a3
3809: if (xot .lt. head) then
3810: uexact = 0.0d+0
3811: elseif (xot .gt. sspd) then
3812: uexact = 0.0d+0
3813: elseif (xot .gt. tail) then
3814: uexact = v2
3815: else
3816: uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
3817: endif
3819: endif
3822: if (debug) then
3823: write(*,*)
3824: !VAM write(*,*) 'x = ',x,' t = ',t
3825: write(*,*) 'uexact = ',uexact
3826: write(*,*)
3827: endif
3829: return
3830: end
3831: double precision function upwind(fw, fe, vp)
3832: !23456789012345678901234567890123456789012345678901234567890123456789012
3833: !
3834: ! function upwind
3835: !
3836: ! this function computes the upwind face value
3837: !
3838: !23456789012345678901234567890123456789012345678901234567890123456789012
3841: !#######################################################################
3843: implicit none
3845: double precision fw, fe, vp
3847: !#######################################################################
3849: if (vp .gt. 0.0) then
3850: upwind = vp * fw
3851: else
3852: upwind = vp * fe
3853: endif
3855: return
3856: end
3857: subroutine uval2
3859: implicit none
3861: double precision prat, grat1, grat2, arat, xnum
3864: logical debug
3866: include 'ex74ftube.h'
3868: debug = .false.
3870: prat = p2/p1
3871: grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
3872: grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
3873: arat = a1/gamma
3875: xnum = sqrt ( grat2 / (prat + grat1) )
3877: v2 = arat * (prat - 1.0d+0) * xnum
3879: if (debug) then
3880: write(*,*)
3881: write(*,*) 'v2 = ',v2
3882: endif
3884: return
3885: end
3886: subroutine val3
3888: implicit none
3890: double precision prat, rpow, epow, p3t
3893: logical debug
3895: include 'ex74ftube.h'
3897: debug = .false.
3900: p3 = p2
3902: prat = p3/p4
3904: rpow = 1.0d+0 / gamma
3906: r3 = r4 * ( prat ** rpow )
3908: epow = (gamma - 1.0d+0) / gamma
3910: e3 = e4 * ( (p3/p4) ** epow )
3912: p3t = (gamma - 1.0d+0) * r3 * e3
3914: a3 = sqrt(gamma*p3/r3)
3916: if (debug) then
3917: write(*,*)
3918: write(*,*) 'a3 = ',a3
3919: write(*,*) 'r3 = ',r3
3920: write(*,*) 'e3 = ',e3
3921: write(*,*) 'p3 = ',p3
3922: write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
3923: write(*,*)
3924: endif
3926: return
3927: end
3928: subroutine wval
3930: implicit none
3932: double precision prat, grat, xnum
3935: logical debug
3937: include 'ex74ftube.h'
3939: debug = .false.
3941: prat = p2/p1
3942: grat = (gamma + 1.0d+0) / (2.0d+0 * gamma)
3944: xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
3945:
3946: sspd = a1 * sqrt(xnum)
3947:
3950: if (debug) then
3951: write(*,*)
3952: write(*,*) 'sspd = ',sspd
3953: endif
3955: return
3956: end