Actual source code: ex74f.F90
petsc-3.7.5 2017-01-01
1: program radhyd
2: ! "$Id: ex4f.F,v 1.39 1999/03/10 19:29:25 Vince Mousseau $";
3: !
4: ! Description: This example solves a nonlinear system on 1 processor with SNES.
5: ! We solve the Euler equations in one dimension.
6: ! The command line options include:
7: ! -dt <dt>, where <dt> indicates time step
8: ! -mx <xg>, where <xg> = number of grid points in the x-direction
9: ! -nstep <nstep>, where <nstep> = number of time steps
10: ! -debug <ndb>, where <ndb> = 0) no debug 1) debug
11: ! -pcnew <npc>, where <npc> = 0) no preconditioning 1) rad preconditioning
12: ! -probnum <probnum>, where <probnum> = 1) cyclic Riesner 2) dam break
13: ! -ihod <ihod>, where <ihod> = 1) upwind 2) quick 3) godunov
14: ! -ientro <ientro>, where <ientro> = 0) basic 1) entropy fix 2) hlle
15: ! -theta <theta>, where <theta> = 0-1 0-explicit 1-implicit
16: ! -hnaught <hnaught>, where <hnaught> = height of left side
17: ! -hlo <hlo>, where <hlo> = hieght of right side
18: ! -ngraph <ngraph>, where <ngraph> = number of time steps between graphics
19: ! -damfac <damfac>, where <damfac> = fractional downward change in hight
20: ! -dampit <ndamp>, where <ndamp> = 1 turn damping on
21: ! -gorder <gorder>, where <gorder> = spatial oerder of godunov
22: !
23: !
24: ! --------------------------------------------------------------------------
25: !
26: ! Shock tube example
27: !
28: ! In this example the application context is a Fortran integer array:
29: ! ctx(1) = shell preconditioner pressure matrix contex
30: ! ctx(2) = semi implicit pressure matrix
31: ! ctx(4) = xold - old time values need for time advancement
32: ! ctx(5) = mx - number of control volumes
33: ! ctx(6) = N - total number of unknowns
34: !
35: ! --------------------------------------------------------------------------
37: implicit none
39: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: ! Include files
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: !
43: #include petsc/finclude/petscdef.h
44: #include petsc/finclude/petsc.h
46: #include "ex74fcomd.h"
47: #include "ex74ftube.h"
49: !
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Variable declarations
52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: !
54: ! Variables:
55: ! snes - nonlinear solver
56: ! x,r - solution, residual vectors
57: ! J - Jacobian matrix
58: ! its - iterations for convergence
59: ! dt - time step size
60: ! draw - drawing context
61: !
62: PetscFortranAddr ctx(6)
63: integer rank,size,nx,ny
64: SNES snes
65: KSP ksp
66: PC pc
67: Vec x,r
68: PetscViewer view0, view1, view2, view3, view4
69: Mat Psemi
70: integer flg, N, ierr, ngraph
71: integer nstep, ndt, i
72: integer its, lits, totits, totlits
73: integer ndb, npc, ndamp, nwilson, ndtcon
74: double precision plotim
76: double precision krtol,katol,kdtol
77: double precision natol,nrtol,nstol
78: integer kmit,nmit,nmf
81: ! Note: Any user-defined Fortran routines (such as FormJacobian)
82: ! MUST be declared as external.
84: external FormFunction, FormInitialGuess,FormDt,PCRadSetUp, PCRadApply, FormGraph,FormDampit
85: double precision eos
87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: ! Initialize program
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: open (unit=87,file='Dt.out',status='unknown')
94: ! start PETSc
95: !
96: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
97: call PetscOptionsSetValue(PETSC_NULL_OBJECT,'-snes_mf','true',ierr)
98: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
99: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
101: if (size .ne. 1) then
102: if (rank .eq. 0) then
103: write(6,*) 'This is a uniprocessor example only!'
104: endif
105: SETERRQ(PETSC_COMM_WORLD,1,0,' ')
106: endif
108: ! Initialize problem parameters
110: debug = .false.
111: dampit = .false.
112: wilson = .true.
113: dtcon = .true.
114: pcnew = .true.
115: dtmax = 1.0d+2
116: dtmin = 1.00d-12
117: dt = 1.0d-2
118: mx = 100
119: nstep = 50
120: probnum = 1
121: gorder = 1
123: tfinal = 1.0d+0
124: tplot = 0.2d+0
125: dtgrow = 1.05d+0
126: tcscal = 0.5d+0
127: hcscal = 0.5d+0
129: ihod = 3
130: ientro = 1
131: theta = 1.00d+0
132: pi = 3.14159d+0
134: zero = 0.0
135: ngraph = 10
137: ndb = 0
138: npc = 1
140: damfac = 0.9d+0
142: gamma = 1.25d+0
143: csubv = 1.0d+0 / (gamma - 1.0d+0)
145: v1 = 0.0d+0
146: v4 = 0.0d+0
148: e1 = 1.0d+0
149: e4 = 1.0d+0
151: r1 = 1.0d+0
152: r4 = 2.0d+0
154: ru1 = r1 * v1
155: ru4 = r4 * v4
157: et1 = r1 * ( (0.5d+0 * v1 * v1) + e1 )
158: et4 = r4 * ( (0.5d+0 * v4 * v4) + e4 )
160: p1 = eos(r1,ru1,et1)
161: p4 = eos(r4,ru4,et4)
163: a1 = sqrt(gamma*p1/r1)
164: a4 = sqrt(gamma*p4/r4)
166: erg0 = 1.0d+2
167: kappa0 = 1.0d+0
168: kappaa = -2.0d+0
169: kappab = 13.0d+0 / 2.0d+0
171: !
172: ! load the command line options
173: !
174: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
175: & '-dt',dt,flg,ierr)
176: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
177: & '-mx',mx,flg,ierr)
178: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
179: & '-nstep',nstep,flg,ierr)
180: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
181: & '-debug',ndb,flg,ierr)
182: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
183: & '-pcnew',npc,flg,ierr)
184: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
185: & '-ihod',ihod,flg,ierr)
186: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
187: & '-ientro',ientro,flg,ierr)
188: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
189: & '-theta',theta,flg,ierr)
190: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
191: & '-ngraph',ngraph,flg,ierr)
192: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
193: & '-damfac',damfac,flg,ierr)
194: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
195: & '-dampit',ndamp,flg,ierr)
196: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
197: & '-wilson',nwilson,flg,ierr)
198: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
199: & '-gorder',gorder,flg,ierr)
200: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
201: & '-probnum',probnum,flg,ierr)
202: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
203: & '-kappa0',kappa0,flg,ierr)
204: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
205: & '-erg0',erg0,flg,ierr)
206: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
207: & '-dtcon',ndtcon,flg,ierr)
208: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
209: & '-tfinal',tfinal,flg,ierr)
210: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
211: & '-tplot',tplot,flg,ierr)
212: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
213: & '-dtgrow',dtgrow,flg,ierr)
214: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
215: & '-tcscal',tcscal,flg,ierr)
216: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
217: & '-hcscal',hcscal,flg,ierr)
218: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
219: & '-dtmax',dtmax,flg,ierr)
221: if (ndamp .eq. 1) then
222: dampit = .true.
223: endif
225: if (nwilson .eq. 0) then
226: wilson = .false.
227: endif
229: if (ndb .eq. 1) then
230: debug = .true.
231: endif
233: if (npc .eq. 0) then
234: pcnew = .false.
235: endif
237: if (ndtcon .eq. 0) then
238: dtcon = .false.
239: endif
241: !CVAM if (dt .ge. dtmax .or. dt .le. dtmin) then
242: !CVAM if (rank .eq. 0) write(6,*) 'DT is out of range'
243: !CVAM SETERRA(1,0,' ')
244: !CVAM endif
246: N = mx*neq
248: ctx(5) = mx
249: ctx(6) = N
251: if (debug) then
252: write(*,*) 'mx = ',mx
253: endif
257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: ! Create nonlinear solver context
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
261: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: ! Create vector data structures; set function evaluation routine
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267: call MatCreateSeqAIJ(PETSC_COMM_WORLD,mx,mx,10,PETSC_NULL_INTEGER,ctx(2),ierr)
269: if (debug) then
270: call MatGetSize(ctx(2),nx,ny,ierr)
271: write(*,*) 'number of rows = ',nx,' number of col = ',ny
272: endif
273: !
274: ! full size vectors
275: !
276: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,x,ierr)
277: call VecSetFromOptions(x,ierr)
278: call VecDuplicate(x,r,ierr)
279: call VecDuplicate(x,ctx(4),ierr)
280: !
281: ! set grid
282: !
283: dx = 2.0d+0/dfloat(mx)
284: xl0 = -1.0d+0 -(0.5d+0 * dx)
286: if (debug) then
287: write(*,*) 'dx = ',dx
288: endif
289:
291: ! Set function evaluation routine and vector. Whenever the nonlinear
292: ! solver needs to evaluate the nonlinear function, it will call this
293: ! routine.
294: ! - Note that the final routine argument is the user-defined
295: ! context that provides application-specific data for the
296: ! function evaluation routine.
298: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
300: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301: ! Customize nonlinear solver; set runtime options
302: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
306: call SNESSetFromOptions(snes,ierr)
307: !
308: ! set the line search function to damp the newton update.
309: !
310: if (dampit) then
311: call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
312: endif
313: !
314: ! get the linear solver info from the nonlinear solver
315: !
317: call SNESGetKSP(snes,ksp,ierr)
318: call KSPGetPC(ksp,pc,ierr)
320: call KSPGetTolerances(ksp,krtol,katol,kdtol,kmit,ierr)
321: call SNESGetTolerances(snes,natol,nrtol,nstol,nmit,nmf,ierr)
323: write(*,*)
324: write(*,*)
325: write(*,*) 'Linear solver'
326: write(*,*)
327: write(*,*) 'rtol = ',krtol
328: write(*,*) 'atol = ',katol
329: write(*,*) 'dtol = ',kdtol
330: write(*,*) 'maxits = ',kmit
331: write(*,*)
332: write(*,*)
333: write(*,*) 'Nonlinear solver'
334: write(*,*)
335: write(*,*) 'rtol = ',nrtol
336: write(*,*) 'atol = ',natol
337: write(*,*) 'stol = ',nstol
338: write(*,*) 'maxits = ',nmit
339: write(*,*) 'max func = ',nmf
340: write(*,*)
341: write(*,*)
343: !
344: ! Build shell based preconditioner if flag set
345: !
346: if (pcnew) then
347: call PCSetType(pc,PCSHELL,ierr)
348: call PCShellSetContext(pc,ctx,ierr)
349: call PCShellSetSetUpCtx(pc,PCRadSetUp,ierr)
350: call PCShellSetApplyCtx(pc,PCRadApply,ctx,ierr)
351: endif
353: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
355: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: ! Evaluate initial guess; then solve nonlinear system.
357: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358: !
359: ! initial counters
360: !
361: time = 0.0d+0
362: plotim = 0.0d+0
363: totits = 0
364: totlits = 0
366: ! Note: The user should initialize the vector, x, with the initial guess
367: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
368: ! to employ an initial guess of zero, the user should explicitly set
369: ! this vector to zero by calling VecSet().
371: call FormInitialGuess(x,ierr)
372: !
373: ! open a window to plot results
374: !
375: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'density',0,0,300,300,view0,ierr)
376: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'velocity',320,0,300,300,view1,ierr)
377: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'total energy',640,0,300,300,view2,ierr)
378: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'temperature',0,360,300,300,view3,ierr)
379: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'pressure',320,360,300,300,view4,ierr)
380: !
381: ! graph initial conditions
382: !
383: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
384: !
385: ! copy x into xold
386: !
387: call VecCopy(x,ctx(4),ierr)
388: call FormDt(snes,x,ctx,ierr)
389: !
390: !################################
391: !
392: ! TIME STEP LOOP BEGIN
393: !
394: !################################
395: !
396: ndt = 0
398: 10 if ( (ndt .le. nstep) .and. ((time + 1.0d-10) .lt. tfinal) ) then
400: if (debug) then
401: write(*,*)
402: write(*,*) 'start of time loop'
403: write(*,*)
404: write(*,*) 'ndt = ',ndt
405: write(*,*) 'nstep = ',nstep
406: write(*,*) 'time = ',time
407: write(*,*) 'tfinal = ',tfinal
408: write(*,*)
409: endif
411: ndt = ndt + 1
412: !
413: ! increment time
414: !
415: time = time + dt
416: plotim = plotim + dt
417: !
418: ! call the nonlinear solver
419: !
420: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
421: !
422: ! get the number of linear iterations used by the nonlinear solver
423: !
424: call SNESGetLinearSolveIterations(snes,lits,ierr)
425: call SNESGetIterationNumber(snes,its,ierr)
427: if (debug) then
428: write(*,*) 'in radhyd ',ndt,'x'
429: call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
430: endif
431: !
432: ! increment the counters
433: !
434: totits = totits + its
435: totlits = totlits + lits
436: !
437: ! Compute new time step
438: !
439: call FormDt(snes,x,ctx,ierr)
440: !
441: ! Draw contour plot of solution
442: !
443: if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
444:
445: plotim = plotim - tplot
448: if (rank .eq. 0) then
449: write(6,100) totits,totlits,ndt,dt,time
450: endif
451: 100 format('Newt = ',i7,' lin =',i7,' step =',i7,' dt = ',e8.3,' time = ',e10.4)
452: !
453: ! graph state conditions
454: !
455: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
457: endif
458: !
459: ! copy x into xold
460: !
461: call VecCopy(x,ctx(4),ierr)
464: goto 10
466: endif
467: !
468: !################################
469: !
470: ! TIME STEP LOOP END
471: !
472: !################################
473: !
475: !
476: ! graph final conditions
477: !
478: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
481: write(*,*)
482: write(*,*)
483: write(*,*) 'total Newton iterations = ',totits
484: write(*,*) 'total linear iterations = ',totlits
485: write(*,*) 'Average Newton per time step = ', dble(totits)/dble(ndt)
486: write(*,*) 'Average Krylov per Newton = ', dble(totlits)/dble(totits)
487: write(*,*)
488: write(*,*)
490: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
491: ! Free work space. All PETSc objects should be destroyed when they
492: ! are no longer needed.
493: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
496: call MatDestroy(ctx(2),ierr)
497: call VecDestroy(x,ierr)
498: call VecDestroy(ctx(4),ierr)
499: call VecDestroy(r,ierr)
500: call SNESDestroy(snes,ierr)
501: call PetscViewerDestroy(view0,ierr)
502: call PetscViewerDestroy(view1,ierr)
503: call PetscViewerDestroy(view2,ierr)
504: call PetscViewerDestroy(view3,ierr)
505: call PetscViewerDestroy(view4,ierr)
507: call PCDestroy(ctx(1),ierr)
509: call PetscFinalize(ierr)
511: close(87)
513: stop
514: end
515: subroutine ApplicationDampit(x,deltx,w,ierr)
516: ! ---------------------------------------------------------------------
517: !
518: ! ApplicationDampit - Damps the newton update, called by
519: ! the higher level routine FormDampit().
520: !
521: ! Input Parameter:
522: ! x - current iterate
523: ! deltx - update
524: !
525: ! Output Parameters:
526: ! x - new iterate
527: ! ierr - error code
528: !
529: ! Notes:
530: ! This routine only damps the density. May want to add energy
531: ! in the future
532: !
534: implicit none
536: ! Common blocks:
537: #include "ex74fcomd.h"
539: ! Input/output variables:
540: PetscScalar x(mx*neq), deltx(mx*neq), w(mx*neq)
541: integer ierr
543: ! Local variables:
544: double precision facmin, fac, newx, xmin, se, dse
545: double precision u,en,rn,run
546: integer i, jr, jru, je
548: 0
550: if (debug) then
551: write(*,*) 'begin damping'
552: do i = 1,mx*neq
553: write(*,*)'x(',i,') = ',x(i)
554: enddo
555: write(*,*)
556: do i = 1,mx*neq
557: write(*,*)'deltx(',i,') = ',deltx(i)
558: enddo
559: endif
561: facmin = 1.0d+0
562: !
563: ! set the scale factor
564: !
565: do i=1,mx
566: !
567: ! set pointers
568: !
569: jr = (neq*i) - 2
570: jru = (neq*i) - 1
571: je = (neq*i)
572: !
573: ! make sure dencity stayes positive
574: !
575: newx = x(jr) - deltx(jr)
576: xmin = damfac * x(jr)
578: if (newx .lt. xmin) then
579: fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
580: if (fac .lt. facmin) then
581: if (debug) then
582: write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
583: endif
584: facmin = fac
585: endif
586: endif
587: !
588: ! make sure Total energy stayes positive
589: !
590: newx = x(je) - deltx(je)
591: xmin = damfac * x(je)
593: if (newx .lt. xmin) then
594: fac = (1.0d+0 - damfac)*x(je)/deltx(je)
595: if (fac .lt. facmin) then
596: if (debug) then
597: write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
598: endif
599: facmin = fac
600: endif
601: endif
602: !
603: ! make sure specific internal energy stayes positive
604: !
605:
606: u = x(jru)/x(jr)
607: se = (x(je)/x(jr)) - (0.5d+0 * u * u)
609: en = x(je) - deltx(je)
610: rn = x(jr) - deltx(jr)
611: run = x(jru) - deltx(jru)
613: dse = se - ( (en/rn) - (0.5d+0 * (run/rn) * (run/rn)) )
616: newx = se - dse
617: xmin = damfac * se
619: if (newx .lt. xmin) then
620: fac = (1.0d+0 - damfac) * se / dse
621: if (fac .lt. facmin) then
622: if (debug) then
623: write(*,*) 'se',i, damfac,facmin,fac,se,dse
624: endif
625: facmin = fac
626: endif
627: endif
629: enddo
630: !
631: ! write out warning
632: !
633: if (facmin .lt. 1.0d+0) then
634: write(*,*) 'facmin = ',facmin, damfac,time
635: !
636: ! scale the vector
637: !
638: do i=1,neq*mx
639: w(i) = x(i) - (facmin * deltx(i))
640: enddo
641: else
642: do i=1,neq*mx
643: w(i) = x(i) - deltx(i)
644: enddo
645: endif
647: if (debug) then
648: write(*,*) 'end damping'
649: do i = 1,mx*neq
650: write(*,*) 'w(',i,') = ',w(i)
651: enddo
652: endif
654: return
655: end
656: subroutine ApplicationDt(x,xold,ierr)
657: ! ---------------------------------------------------------------------
658: !
659: ! ApplicationDt - compute CFL numbers. Called by
660: ! the higher level routine FormDt().
661: !
662: ! Input Parameter:
663: ! x - local vector data
664: !
665: ! Output Parameters:
666: ! ierr - error code
667: !
668: ! Notes:
669: ! This routine uses standard Fortran-style computations over a 2-dim array.
670: !
672: implicit none
674: ! Common blocks:
675: #include "ex74fcomd.h"
676: #include "ex74ftube.h"
678: ! Input/output variables:
679: PetscScalar x(mx*neq), xold(mx*neq)
680: integer ierr
682: ! Local variables:
683: integer i, jr, jru, je
684: !
685: ! new
686: !
687: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw
688: double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
689: double precision tcflg, hcflg, dtt, dth
690: double precision te, tp, tw
691: double precision ue, up, uw
692: double precision see, sep, sew
693: !
694: ! old
695: !
696: double precision rhooee, rhooe, rhoop, rhoow, rhooww
697: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
698: double precision ergoee, ergoe, ergop, ergow, ergoww
699: double precision veloe, velop, velow
700: double precision uop, seop, top
701: double precision dtold, dttype
702: !
703: ! functions
704: !
705: double precision eos
707: dtold = dt
709: 0
711: if (debug) then
712: write(*,*) 'in start dt'
713: do i = 1,mx*neq
714: write(*,*)'x(',i,') = ',x(i)
715: enddo
716: write(*,*) 'tfinal = ',tfinal
717: write(*,*) 'time = ',time
718: write(*,*) 'dt = ',dt
719: write(*,*) 'dtmax = ',dtmax
720: endif
722: sndp = -1.0d+20
723: vradn = 0.0d+0
724: vradd = 0.0d+0
726: !
727: !################################
728: !
729: ! loop over all cells begin
730: !
731: !################################
732: !
733: do i=1,mx
734: !
735: ! set pointers
736: !
737: jr = (neq*i) - 2
738: jru = (neq*i) - 1
739: je = (neq*i)
740: !
741: !
742: ! set scalars
743: !
744: call Setpbc(i,x, rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww,vele, velp, velw)
745: !
746: ! compute temperatures
747: !
748: uw = rhouw / rhow
749: up = rhoup / rhop
750: ue = rhoue / rhoe
752: see = (erge/rhoe) - (0.5d+0 * ue * ue)
753: sep = (ergp/rhop) - (0.5d+0 * up * up)
754: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
756: te = see / csubv
757: tp = sep / csubv
758: tw = sew / csubv
759: !
760: ! compute old temperature
761: !
763: call Setpbc(i,xold,rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
765: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww,veloe, velop, velow, i)
767: uop = rhouop / rhoop
769: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
771: top = seop / csubv
772: !
773: ! compute thermal cfl
774: !
775: vradn = vradn + (abs(tp - top)/dt)
776: vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
777: !
778: ! compute hydro cfl
779: !
781: pressp = eos(rhop, rhoup, ergp)
782: sndp = max(sndp,sqrt( (gamma*pressp) / rhop ))
784: enddo
785: !
786: !################################
787: !
788: ! loop over all cells end
789: !
790: !################################
791: !
793: vrad = vradn / vradd
795: tcfl = (vrad * dt) / dx
796: hcfl = (sndp * dt) / dx
798: dtt = max(dx/vrad,1.0d-7)
799: dtt = tcscal * dtt
801: dth = hcscal * dx / sndp
803: if (.not. dtcon) then
804: dt = min (dth,dtt,dt*dtgrow)
805: if (dt .lt. dtmin) then
806: dt = dtmin
807: endif
808: if (dt .gt. dtmax) then
809: dt = dtmax
810: endif
811: if ( (time + dt) .gt. tfinal) then
812: dt = tfinal - time
813: endif
815: if (dt .eq. dth) then
816: dttype = 1.0d+0
817: elseif (dt .eq. dtt) then
818: dttype = 2.0d+0
819: elseif (dt .eq. dtold*dtgrow) then
820: dttype = 3.0d+0
821: elseif (dt .eq. dtmax) then
822: dttype = 4.0d+0
823: elseif (dt .eq. dtmin) then
824: dttype = 5.0d+0
825: elseif (dt .eq. tfinal - time) then
826: dttype = 6.0
827: else
828: dttype = -1.0d+0
829: endif
831: endif
832:
833:
834: write(87,1000) time,dt,dth/hcscal,dtt/tcscal
835: write(88,1000) time,dttype
837: 1000 format(4(2x,e18.9))
839: if (debug) then
840: write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
841: write(*,*) 'dtt = ',dtt,' dth = ',dth
842: write(*,*) 'tfinal = ',tfinal
843: write(*,*) 'time = ',time
844: write(*,*) 'dt = ',dt
845: write(*,*) 'dtmax = ',dtmax
846: write(*,*)
847: endif
849: return
850: end
851: subroutine ApplicationExact(x,ierr)
852: ! ---------------------------------------------------------------------
853: !
854: ! ApplicationExact - Computes exact solution, called by
855: ! the higher level routine FormExact().
856: !
857: ! Input Parameter:
858: ! x - local vector data
859: !
860: ! Output Parameters:
861: ! x - initial conditions
862: ! ierr - error code
863: !
864: ! Notes:
865: ! This routine uses standard Fortran-style computations over a 1-dim array.
866: !
868: implicit none
870: ! Common blocks:
872: #include "ex74fcomd.h"
874: ! Input/output variables:
875: PetscScalar x(mx)
876: integer ierr
878: ! Local variables:
879: integer i
880: double precision xloc
881: PetscScalar rexact
884: ! Set parameters
886: 0
888: do i = 1,mx
890: xloc = xl0 + (dble(i) * dx)
891: x(i) = rexact(xloc,time)
893: enddo
895: return
896: end
897: subroutine ApplicationFunction(x,f,xold,ierr)
898: ! ---------------------------------------------------------------------
899: !
900: ! ApplicationFunction - Computes nonlinear function, called by
901: ! the higher level routine FormFunction().
902: !
903: ! Input Parameter:
904: ! x - local vector data
905: !
906: ! Output Parameters:
907: ! f - local vector data, f(x)
908: ! ierr - error code
909: !
910: ! Notes:
911: ! This routine uses standard Fortran-style computations over a 2-dim array.
912: !
914: implicit none
916: ! Common blocks:
917: #include "ex74fcomd.h"
919: ! Input/output variables:
920: PetscScalar x(mx*neq), f(mx*neq), xold(mx*neq)
921: integer ierr
923: ! Local variables:
924: integer i, jr, jru, je
925: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw
927: double precision cont, energy, mom
929: 0
931: if (debug) then
932: write(*,*) 'in function'
933: do i = 1,mx*neq
934: write(*,*)'x(',i,') = ',x(i)
935: enddo
936: endif
937: !
938: !################################
939: !
940: ! loop over all cells begin
941: !
942: !################################
943: !
944: do i=1,mx
945: !
946: ! set pointers
947: !
948: jr = (neq*i) - 2
949: jru = (neq*i) - 1
950: je = (neq*i)
951: !
952: !
953: ! set scalars
954: !
955: call Setpbc(i,x,rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw)
956: !
957: ! compute functions
958: !
960: f(jr) = cont(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
963: f(jru) = mom(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
966: f(je) = energy(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, i,xold)
968: if (debug) then
969: write(*,*)
970: write(*,*) i,jr,jru,je,'res,r,ru,e'
971: write(*,*) f(jr),f(jru),f(je)
972: write(*,*)
973: endif
975: enddo
976: !
977: !################################
978: !
979: ! loop over all cells end
980: !
981: !################################
982: !
984: if (debug) then
985: write(*,*) 'in function'
986: do i = 1,mx*neq
987: write(*,*) 'f(',i,') = ',f(i)
988: enddo
989: endif
991: return
992: end
993: subroutine ApplicationInitialGuess(x,ierr)
994: ! ---------------------------------------------------------------------
995: !
996: ! ApplicationInitialGuess - Computes initial approximation, called by
997: ! the higher level routine FormInitialGuess().
998: !
999: ! Input Parameter:
1000: ! x - local vector data
1001: !
1002: ! Output Parameters:
1003: ! x - initial conditions
1004: ! ierr - error code
1005: !
1006: ! Notes:
1007: ! This routine uses standard Fortran-style computations over a 1-dim array.
1008: !
1010: implicit none
1012: ! Common blocks:
1014: #include "ex74fcomd.h"
1015: #include "ex74ftube.h"
1017: ! Input/output variables:
1018: PetscScalar x(mx*neq)
1019: integer ierr
1021: ! Local variables:
1022: integer i, j, jr, jru, je
1023: double precision xloc, re, ee, ve
1024: double precision wid, efloor
1025: PetscScalar uexact, rexact, eexact
1028: !VAM efloor = max(1.0d-1,1.0d-3 * erg0)
1029: efloor = 1.0d-1
1030: !VAM wid = max(1.0d-2,dx)
1031: wid = 1.0d-2
1033: ! Set parameters
1035: 0
1037: do i = 1,mx
1039: jr = (neq*i) - 2
1040: jru = (neq*i) - 1
1041: je = (neq*i)
1043: xloc = xl0 + (dble(i) * dx)
1045: if (probnum .eq. 1) then
1046: re = rexact(xloc,zero)
1047: ve = uexact(xloc,zero)
1048: ee = eexact(xloc,zero)
1049: else
1050: re = 1.0d+0
1051: ve = 0.0d+0
1052: ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1053: endif
1055: x(jr) = re
1056: x(jru) = re * ve
1057: x(je) = re * ( (0.5d+0 * ve * ve) + ee )
1059: if (debug) then
1060: write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1061: 100 format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1062: endif
1064: enddo
1066: call exact0
1067: call eval2
1068: call rval2
1069: call wval
1070: call uval2
1071: v3 = v2
1072: call val3
1074: a1 = sqrt(gamma*p1/r1)
1075: a2 = sqrt(gamma*p2/r2)
1076: a3 = sqrt(gamma*p3/r3)
1077: a4 = sqrt(gamma*p4/r4)
1079: write(*,1000) r1,r2,r3,r4
1080: write(*,2000) p1,p2,p3,p4
1081: write(*,3000) e1,e2,e3,e4
1082: write(*,4000) a1,a2,a3,a4
1083: write(*,*)
1085: 1000 format ('rhos ',4(f12.6))
1086: 2000 format ('pressures ',4(f12.6))
1087: 3000 format ('energies ',4(f12.6))
1088: 4000 format ('sound ',4(f12.6))
1091: return
1092: end
1093: subroutine ApplicationXmgr(x,ivar,ierr)
1094: ! ---------------------------------------------------------------------
1095: !
1096: ! ApplicationXmgr - Sets the Xmgr output called from
1097: ! the higher level routine FormXmgr().
1098: !
1099: ! Input Parameter:
1100: ! x - local vector data
1101: !
1102: ! Output Parameters:
1103: ! x - initial conditions
1104: ! ierr - error code
1105: !
1106: ! Notes:
1107: ! This routine uses standard Fortran-style computations over a 1-dim array.
1108: !
1110: implicit none
1112: ! Common blocks:
1114: #include "ex74fcomd.h"
1116: ! Input/output variables:
1117: PetscScalar x(mx)
1118: integer ivar,ierr
1120: ! Local variables:
1121: integer i
1122: double precision xloc, sum
1123: PetscScalar rexact
1124: integer iplotnum(5)
1125: save iplotnum
1126: character*8 grfile
1129: data iplotnum / -1,-1,-1,-1,-1 /
1133: ! Set parameters
1135: iplotnum(ivar) = iplotnum(ivar) + 1
1136: 0
1138: if (ivar .eq. 1) then
1139: write(grfile,4000) iplotnum(ivar)
1140: 4000 format('Xmgrr',i3.3)
1141: elseif (ivar .eq. 2) then
1142: write(grfile,5000) iplotnum(ivar)
1143: 5000 format('Xmgru',i3.3)
1144: elseif (ivar .eq. 3) then
1145: write(grfile,6000) iplotnum(ivar)
1146: 6000 format('Xmgre',i3.3)
1147: elseif (ivar .eq. 4) then
1148: write(grfile,7000) iplotnum(ivar)
1149: 7000 format('Xmgrt',i3.3)
1150: else
1151: write(grfile,8000) iplotnum(ivar)
1152: 8000 format('Xmgrp',i3.3)
1153: endif
1155: open (unit=44,file=grfile,status='unknown')
1157: do i = 1,mx
1159: xloc = xl0 + (dble(i) * dx)
1160: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1161: write(44,1000) xloc, x(i), rexact(xloc,time)
1162: else
1163: write(44,1000) xloc, x(i)
1164: endif
1166: enddo
1168: 1000 format(3(e18.12,2x))
1169: close(44)
1171: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1172: sum = 0.0d+0
1173: do i = 1,mx
1174: xloc = xl0 + (dble(i) * dx)
1175: sum = sum + (x(i) - rexact(xloc,time)) ** 2
1176: enddo
1177: sum = sqrt(sum)
1179: write(*,*)
1180: write(*,*) 'l2norm of the density error is',sum
1181: write(*,*)
1182: endif
1185: return
1186: end
1187: subroutine FormDampit(snes,ctx,x,f,g,y,w, fnorm,ynorm,gnorm,flag,ierr)
1188: ! ---------------------------------------------------------------------
1189: !
1190: ! FormDampit - damps the Newton update
1191: !
1192: ! Input Parameters:
1193: ! snes - the SNES context
1194: ! x - current iterate
1195: ! f - residual evaluated at x
1196: ! y - search direction (containes new iterate on output)
1197: ! w - work vector
1198: ! fnorm - 2-norm of f
1199: !
1200: ! In this example the application context is a Fortran integer array:
1201: ! ctx(1) = shell preconditioner pressure matrix contex
1202: ! ctx(2) = semi implicit pressure matrix
1203: ! ctx(4) = xold - old time values need for time advancement
1204: ! ctx(5) = mx - number of control volumes
1205: ! ctx(6) = N - total number of unknowns
1206: !
1207: ! Output Parameter:
1208: ! g - residual evaluated at new iterate y
1209: ! y - new iterate (contains search direction on input
1210: ! gnorm - 2-norm of g
1211: ! ynorm - 2-norm of search length
1212: ! flag - set to 0 if the line search succeeds; -1 on failure
1213: !
1214: ! Notes:
1215: ! This routine serves as a wrapper for the lower-level routine
1216: ! "ApplicationDampit", where the actual computations are
1217: ! done using the standard Fortran style of treating the local
1218: ! vector data as a multidimensional array over the local mesh.
1219: ! This routine merely accesses the local vector data via
1220: ! VecGetArray() and VecRestoreArray().
1221: !
1222: implicit none
1224: #include petsc/finclude/petsc.h
1226: ! Input/output variables:
1227: SNES snes
1228: Vec x, f, g, y, w
1229: PetscFortranAddr ctx(*)
1230: PetscScalar fnorm, ynorm, gnorm
1231: integer ierr, flag
1233: ! Common blocks:
1235: #include "ex74fcomd.h"
1237: ! Local variables:
1239: ! Declarations for use with local arrays:
1240: PetscScalar lx_v(0:1), ly_v(0:1), lw_v(0:1)
1241: PetscOffset lx_i, ly_i , lw_i
1243: !
1244: ! set ynorm
1245: !
1246: call VecNorm(y,NORM_2,ynorm,ierr)
1247: !
1248: ! copy x to w
1249: !
1250: call VecCopy(x,w,ierr)
1251: !
1252: ! get pointers to x, y, w
1253: !
1255: call VecGetArray(x,lx_v,lx_i,ierr)
1256: call VecGetArray(y,ly_v,ly_i,ierr)
1257: call VecGetArray(w,lw_v,lw_i,ierr)
1258: !
1259: ! Compute Damping
1260: !
1261: call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1262: !
1263: ! Restore vectors x, y, w
1264: !
1265: call VecRestoreArray(x,lx_v,lx_i,ierr)
1266: call VecRestoreArray(y,ly_v,ly_i,ierr)
1267: call VecRestoreArray(w,lw_v,lw_i,ierr)
1268: !
1269: ! copy w to y
1270: !
1271: call VecCopy(w,y,ierr)
1272: !
1273: ! compute new residual
1274: !
1275: call SNESComputeFunction(snes,y,g,ierr)
1276: call VecNorm(g,NORM_2,gnorm,ierr)
1277: flag = 0
1279: if (debug) then
1280: write(*,*) 'form damp ynorm = ',ynorm
1281: write(*,*) 'gnorm = ',gnorm
1282: endif
1284: return
1285: end
1286: subroutine FormDt(snes,x,ctx,ierr)
1287: ! ---------------------------------------------------------------------
1288: !
1289: ! FormDt - Compute CFL numbers
1290: !
1291: ! Input Parameters:
1292: ! snes - the SNES context
1293: ! x - input vector
1294: !
1295: ! In this example the application context is a Fortran integer array:
1296: ! ctx(1) = shell preconditioner pressure matrix contex
1297: ! ctx(2) = semi implicit pressure matrix
1298: ! ctx(4) = xold - old time values need for time advancement
1299: ! ctx(5) = mx - number of control volumes
1300: ! ctx(6) = N - total number of unknowns
1301: !
1302: !
1303: ! Notes:
1304: ! This routine serves as a wrapper for the lower-level routine
1305: ! "ApplicationDt", where the actual computations are
1306: ! done using the standard Fortran style of treating the local
1307: ! vector data as a multidimensional array over the local mesh.
1308: ! This routine merely accesses the local vector data via
1309: ! VecGetArray() and VecRestoreArray().
1310: !
1311: implicit none
1313: #include petsc/finclude/petsc.h
1315: ! Input/output variables:
1316: SNES snes
1317: Vec x
1318: PetscFortranAddr ctx(*)
1319: integer ierr
1321: ! Common blocks:
1323: #include "ex74fcomd.h"
1325: ! Local variables:
1327: ! Declarations for use with local arrays:
1328: PetscScalar lx_v(0:1)
1329: PetscOffset lx_i
1330: PetscScalar lxold_v(0:1)
1331: PetscOffset lxold_i
1333: !
1334: ! get pointers to x, xold
1335: !
1337: call VecGetArray(x,lx_v,lx_i,ierr)
1338: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1339: !
1340: ! Compute function
1341: !
1342: call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1343: !
1344: ! Restore vectors x, xold
1345: !
1346: call VecRestoreArray(x,lx_v,lx_i,ierr)
1347: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1349: return
1350: end
1351: subroutine FormExact(x,ierr)
1352: ! ---------------------------------------------------------------------
1353: !
1354: ! FormExact - Forms exact solution
1355: !
1356: ! Input Parameter:
1357: ! x - vector
1358: !
1359: ! Output Parameters:
1360: ! x - vector
1361: ! ierr - error code
1362: !
1363: ! Notes:
1364: ! This routine serves as a wrapper for the lower-level routine
1365: ! "ApplicationExact", where the actual computations are
1366: ! done using the standard Fortran style of treating the local
1367: ! vector data as a multidimensional array over the local mesh.
1368: ! This routine merely accesses the local vector data via
1369: ! VecGetArray() and VecRestoreArray().
1370: !
1371: implicit none
1373: #include petsc/finclude/petsc.h
1375: ! Input/output variables:
1376: Vec x
1377: integer ierr
1379: ! Declarations for use with local arrays:
1380: PetscScalar lx_v(0:1)
1381: PetscOffset lx_i
1383: 0
1385: !
1386: ! get a pointer to x
1387: !
1388: call VecGetArray(x,lx_v,lx_i,ierr)
1389: !
1390: ! Compute initial guess
1391: !
1392: call ApplicationExact(lx_v(lx_i),ierr)
1393: !
1394: ! Restore vector x
1395: !
1396: call VecRestoreArray(x,lx_v,lx_i,ierr)
1398: return
1399: end
1400: subroutine FormFunction(snes,x,f,ctx,ierr)
1401: ! ---------------------------------------------------------------------
1402: !
1403: ! FormFunction - Evaluates nonlinear function, f(x).
1404: !
1405: ! Input Parameters:
1406: ! snes - the SNES context
1407: ! x - input vector
1408: !
1409: ! In this example the application context is a Fortran integer array:
1410: ! ctx(1) = shell preconditioner pressure matrix contex
1411: ! ctx(2) = semi implicit pressure matrix
1412: ! ctx(4) = xold - old time values need for time advancement
1413: ! ctx(5) = mx - number of control volumes
1414: ! ctx(6) = N - total number of unknowns
1415: !
1416: ! Output Parameter:
1417: ! f - vector with newly computed function
1418: !
1419: ! Notes:
1420: ! This routine serves as a wrapper for the lower-level routine
1421: ! "ApplicationFunction", where the actual computations are
1422: ! done using the standard Fortran style of treating the local
1423: ! vector data as a multidimensional array over the local mesh.
1424: ! This routine merely accesses the local vector data via
1425: ! VecGetArray() and VecRestoreArray().
1426: !
1427: implicit none
1429: #include petsc/finclude/petsc.h
1431: ! Input/output variables:
1432: SNES snes
1433: Vec x, f
1434: PetscFortranAddr ctx(*)
1435: integer ierr
1437: ! Common blocks:
1439: #include "ex74fcomd.h"
1441: ! Local variables:
1443: ! Declarations for use with local arrays:
1444: PetscScalar lx_v(0:1), lf_v(0:1)
1445: PetscOffset lx_i, lf_i
1446: PetscScalar lxold_v(0:1)
1447: PetscOffset lxold_i
1449: !
1450: ! get pointers to x, f, xold
1451: !
1453: call VecGetArray(x,lx_v,lx_i,ierr)
1454: call VecGetArray(f,lf_v,lf_i,ierr)
1455: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1456: !
1457: ! Compute function
1458: !
1459: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i), lxold_v(lxold_i),ierr)
1460: !
1461: ! Restore vectors x, f, xold
1462: !
1463: call VecRestoreArray(x,lx_v,lx_i,ierr)
1464: call VecRestoreArray(f,lf_v,lf_i,ierr)
1465: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1466: !
1467: ! something to do with profiling
1468: !
1469: call PetscLogFlops(110.d0*mx,ierr)
1471: return
1472: end
1473: subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1474: ! ---------------------------------------------------------------------
1475: !
1476: ! FormGraph - Forms Graphics output
1477: !
1478: ! Input Parameter:
1479: ! x - vector
1480: ! time - scalar
1481: !
1482: ! Output Parameters:
1483: ! ierr - error code
1484: !
1485: ! Notes:
1486: ! This routine serves as a wrapper for the lower-level routine
1487: ! "ApplicationXmgr", where the actual computations are
1488: ! done using the standard Fortran style of treating the local
1489: ! vector data as a multidimensional array over the local mesh.
1490: ! This routine merely accesses the local vector data via
1491: ! VecGetArray() and VecRestoreArray().
1492: !
1493: implicit none
1495: #include petsc/finclude/petsc.h
1497: #include "ex74fcomd.h"
1498: #include "ex74ftube.h"
1500: ! Input/output variables:
1501: Vec x
1502: integer ierr
1504: ! Declarations for use with local arrays:
1505: IS rfrom, rto, rufrom, ruto, efrom, eto
1506: Vec rval
1507: Vec uval
1508: Vec ruval
1509: Vec eval
1510: Vec seval
1511: Vec pval
1512: Vec kval
1513: Vec tval
1514: Vec steval
1515: VecScatter scatter
1516: PetscViewer view0, view1, view2, view3, view4
1517: double precision minus1, l2err, gm1, csubvi
1520: csubvi = 1.0d+0 / csubv
1521: gm1 = gamma - 1.0d+0
1522: 0
1523: minus1 = -1.0d+0
1524: !
1525: ! graphics vectors
1526: !
1527: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,rval,ierr)
1528: call VecSetFromOptions(rval,ierr)
1529: call VecDuplicate(rval,uval,ierr)
1530: call VecDuplicate(rval,ruval,ierr)
1531: call VecDuplicate(rval,eval,ierr)
1532: call VecDuplicate(rval,seval,ierr)
1533: call VecDuplicate(rval,pval,ierr)
1534: call VecDuplicate(rval,kval,ierr)
1535: call VecDuplicate(rval,tval,ierr)
1536: call VecDuplicate(rval,steval,ierr)
1537: !
1538: ! create index sets for vector scatters
1539: !
1540: call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1541: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, rto, ierr)
1542: call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1543: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, ruto, ierr)
1544: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1545: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, eto, ierr)
1547: !
1548: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1549: !
1550: !
1551: ! load rval from x
1552: !
1553: call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1554: call VecScatterBegin(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1555: call VecScatterEnd(scatter,x,rval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1556: call VecScatterDestroy(scatter,ierr)
1557: !
1558: ! plot rval vector
1559: !
1560: call VecView(rval,view0,ierr)
1561: !
1562: ! make xmgr plot of rval
1563: !
1564: call FormXmgr(rval,1,ierr)
1565: !
1566: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1567: !
1568: !
1569: ! load eval from x
1570: !
1571: call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1572: call VecScatterBegin(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1573: call VecScatterEnd(scatter,x,eval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1574: call VecScatterDestroy(scatter,ierr)
1575: !
1576: ! plot eval vector
1577: !
1578: call VecView(eval,view2,ierr)
1579: !
1580: ! make xmgr plot of eval
1581: !
1582: call FormXmgr(eval,3,ierr)
1583: !
1584: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1585: !
1587: !
1588: ! load ruval from x
1589: !
1590: call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1591: call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,SCATTER_FORWARD,ierr)
1592: call VecScatterEnd(scatter,x,ruval,INSERT_VALUES, SCATTER_FORWARD,ierr)
1593: call VecScatterDestroy(scatter,ierr)
1594: !
1595: ! create u = ru / r
1596: !
1597: call VecPointwiseDivide(uval,ruval,rval,ierr)
1598: !
1599: ! plot uval vector
1600: !
1601: call VecView(uval,view1,ierr)
1602: !
1603: ! make xmgr plot of uval
1604: !
1605: call FormXmgr(uval,2,ierr)
1607: !
1608: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1609: !
1611: call VecPointwiseMult(kval,uval,uval,ierr)
1612: call VecScale(kval,0.5d+0,ierr)
1614: call VecPointwiseDivide(steval,eval,rval,ierr)
1615: call VecWAXPY(seval,-1.0d+0,kval,steval,ierr)
1617: call VecCopy(seval,tval,ierr)
1618: call VecScale(tval,csubvi,ierr)
1620: !
1621: ! plot tval vector
1622: !
1623: call VecView(tval,view3,ierr)
1624: !
1625: ! make xmgr plot of tval
1626: !
1627: call FormXmgr(tval,4,ierr)
1629: !
1630: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1631: !
1633: call VecPointwiseMult(pval,rval,seval,ierr)
1634: call VecScale(pval,gm1,ierr)
1635: !
1636: ! plot pval vector
1637: !
1638: call VecView(pval,view4,ierr)
1639: !
1640: ! make xmgr plot of pval
1641: !
1642: call FormXmgr(pval,5,ierr)
1643: !
1644: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1645: !
1651: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1652: ! Free work space. All PETSc objects should be destroyed when they
1653: ! are no longer needed.
1654: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1656: call VecDestroy(rval, ierr)
1657: call VecDestroy(uval, ierr)
1658: call VecDestroy(ruval,ierr)
1659: call VecDestroy(eval, ierr)
1660: call VecDestroy(seval, ierr)
1661: call VecDestroy(pval, ierr)
1662: call VecDestroy(kval, ierr)
1663: call VecDestroy(tval, ierr)
1664: call VecDestroy(steval, ierr)
1666: call ISDestroy(rfrom, ierr)
1667: call ISDestroy(rto, ierr)
1669: call ISDestroy(rufrom,ierr)
1670: call ISDestroy(ruto, ierr)
1672: call ISDestroy(efrom, ierr)
1673: call ISDestroy(eto, ierr)
1676: return
1677: end
1678: subroutine FormInitialGuess(x,ierr)
1679: ! ---------------------------------------------------------------------
1680: !
1681: ! FormInitialGuess - Forms initial approximation.
1682: !
1683: ! Input Parameter:
1684: ! x - vector
1685: !
1686: ! Output Parameters:
1687: ! x - vector
1688: ! ierr - error code
1689: !
1690: ! Notes:
1691: ! This routine serves as a wrapper for the lower-level routine
1692: ! "ApplicationInitialGuess", where the actual computations are
1693: ! done using the standard Fortran style of treating the local
1694: ! vector data as a multidimensional array over the local mesh.
1695: ! This routine merely accesses the local vector data via
1696: ! VecGetArray() and VecRestoreArray().
1697: !
1698: implicit none
1700: #include petsc/finclude/petsc.h
1702: ! Input/output variables:
1703: Vec x
1704: integer ierr
1706: ! Declarations for use with local arrays:
1707: PetscScalar lx_v(0:1)
1708: PetscOffset lx_i
1710: 0
1712: !
1713: ! get a pointer to x
1714: !
1715: call VecGetArray(x,lx_v,lx_i,ierr)
1716: !
1717: ! Compute initial guess
1718: !
1719: call ApplicationInitialGuess(lx_v(lx_i),ierr)
1720: !
1721: ! Restore vector x
1722: !
1723: call VecRestoreArray(x,lx_v,lx_i,ierr)
1725: return
1726: end
1727: subroutine FormXmgr(x,ivar,ierr)
1728: ! ---------------------------------------------------------------------
1729: !
1730: ! FormXmgr - Forms Xmgr output
1731: !
1732: ! Input Parameter:
1733: ! x - vector
1734: !
1735: ! Output Parameters:
1736: ! x - vector
1737: ! ierr - error code
1738: !
1739: ! Notes:
1740: ! This routine serves as a wrapper for the lower-level routine
1741: ! "ApplicationXmgr", where the actual computations are
1742: ! done using the standard Fortran style of treating the local
1743: ! vector data as a multidimensional array over the local mesh.
1744: ! This routine merely accesses the local vector data via
1745: ! VecGetArray() and VecRestoreArray().
1746: !
1747: implicit none
1749: #include petsc/finclude/petsc.h
1751: ! Input/output variables:
1752: Vec x
1753: integer ivar,ierr
1755: ! Declarations for use with local arrays:
1756: PetscScalar lx_v(0:1)
1757: PetscOffset lx_i
1759: 0
1761: !
1762: ! get a pointer to x
1763: !
1764: call VecGetArray(x,lx_v,lx_i,ierr)
1765: !
1766: ! make the graph
1767: !
1768: call ApplicationXmgr(lx_v(lx_i),ivar,ierr)
1769: !
1770: ! Restore vector x
1771: !
1772: call VecRestoreArray(x,lx_v,lx_i,ierr)
1774: return
1775: end
1776: subroutine PCRadApply(pc,ctx,x,y,ierr)
1777: ! -------------------------------------------------------------------
1778: !
1779: ! PCRadApply - This routine demonstrates the use of a
1780: ! user-provided preconditioner.
1781: !
1782: ! Input Parameters:
1783: ! dummy - optional user-defined context, not used here
1784: ! x - input vector
1785: ! In this example the shell preconditioner application context
1786: ! is a Fortran integer array:
1787: ! ctx(1) = shell preconditioner pressure matrix contex
1788: ! ctx(2) = semi implicit pressure matrix
1789: ! ctx(4) = xold - old time values need for time advancement
1790: ! ctx(5) = mx - number of control volumes
1791: ! ctx(6) = N - total number of unknowns
1792: !
1793: ! Output Parameters:
1794: ! y - preconditioned vector
1795: ! ierr - error code (nonzero if error has been detected)
1796: !
1797: ! Notes:
1798: ! This code implements the Jacobi preconditioner plus the
1799: ! SOR preconditioner
1800: !
1802: implicit none
1804: #include petsc/finclude/petsc.h
1806: #include "ex74fcomd.h"
1807: !
1808: ! Input
1809: !
1810: PC pc
1811: PetscFortranAddr ctx(*)
1812: Vec x, y
1813: integer ierr
1814: !
1815: ! Local
1816: !
1817: IS defrom, deto
1818: Vec de, rese
1819: VecScatter scatter
1820: PetscScalar lde_v(0:1),lrese_v(0:1)
1821: PetscOffset lde_i, lrese_i
1822: !
1823: ! Identity preconditioner
1824: !
1825: call VecCopy(x,y,ierr)
1826: !
1827: ! if kappa0 not equal to zero then precondition the radiation diffusion
1828: !
1829: if (kappa0 .ne. 0.0d+0) then
1830:
1832: !
1833: ! Create needed vectors
1834: !
1835: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,de,ierr)
1836: call VecSetFromOptions(de,ierr)
1837: call VecDuplicate(de,rese,ierr)
1838: !
1839: ! create index sets for scatters
1840: !
1841: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,defrom,ierr)
1842: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1,deto,ierr)
1843: !
1844: ! load rese from x
1845: !
1846: call VecScatterCreate(x,defrom,rese,deto,scatter,ierr)
1847: call VecScatterBegin(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1848: call VecScatterEnd(scatter,x,rese,INSERT_VALUES, SCATTER_FORWARD,ierr)
1849: call VecScatterDestroy(scatter,ierr)
1850: !
1851: ! apply preconditioner
1852: !
1853: call PCApply(ctx(1),rese,de,ierr)
1855: if (debug) then
1856: write(*,*) 'PCRadApply dh is'
1857: call VecView(de,PETSC_VIEWER_STDOUT_SELF,ierr)
1858: endif
1859: !
1860: ! load de into y
1861: !
1862: call VecScatterCreate(de,deto,y,defrom,scatter,ierr)
1863: call VecScatterBegin(scatter,de,y,INSERT_VALUES, SCATTER_FORWARD,ierr)
1864: call VecScatterEnd(scatter,de,y,INSERT_VALUES,SCATTER_FORWARD,ierr)
1865: call VecScatterDestroy(scatter,ierr)
1867: if (debug) then
1868: write(*,*) 'PCRadApply y is'
1869: call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr)
1870: endif
1872: call VecDestroy(de,ierr)
1873: call VecDestroy(rese,ierr)
1875: call ISDestroy(defrom,ierr)
1876: call ISDestroy(deto,ierr)
1878: endif
1881: return
1882: end
1883: subroutine PCRadSetUp(pc,ctx,ierr)
1884: !
1885: ! PCRadSetUp - This routine sets up a user-defined
1886: ! preconditioner context.
1887: !
1888: ! Input Parameters:
1889: ! In this example the shell preconditioner application context
1890: ! is a Fortran integer array:
1891: ! ctx(1) = shell preconditioner pressure matrix contex
1892: ! ctx(2) = semi implicit pressure matrix
1893: ! ctx(4) = xold - old time values need for time advancement
1894: ! ctx(5) = mx - number of control volumes
1895: ! ctx(6) = N - total number of unknowns
1896: !
1897: ! Output Parameter:
1898: ! ierr - error code (nonzero if error has been detected)
1899: !
1900: ! Notes:
1901: ! In this example, we define the shell preconditioner to be Jacobi
1902: ! method. Thus, here we create a work vector for storing the reciprocal
1903: ! of the diagonal of the preconditioner matrix; this vector is then
1904: ! used within the routine PCRadApply().
1905: !
1907: implicit none
1909: #include petsc/finclude/petsc.h
1911: #include "ex74fcomd.h"
1912: !
1913: ! Input
1914: !
1915: PC pc
1916: PetscFortranAddr ctx(*)
1917: integer ierr
1918: !
1919: ! Local
1920: !
1921: Vec eold
1922:
1923: PetscScalar le_v(0:1)
1924: PetscOffset le_i
1925:
1926: !
1927: ! create vector
1928: !
1929: call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,mx,eold,ierr)
1930: call VecSetFromOptions(eold,ierr)
1931: !
1932: ! set up the matrix based on xold
1933: !
1934: call Setmat(ctx,ierr)
1935: !
1936: ! set up the preconditioner
1937: !
1938: call PCDestroy(ctx(1),ierr)
1939: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
1940: !VAM call PCSetType(ctx(1),PCJACOBI,ierr)
1941: call PCSetType(ctx(1),PCLU,ierr)
1942: ! call PCSetVector(ctx(1),eold,ierr)
1943: call PCSetOperators(ctx(1),ctx(2),ctx(2), ierr)
1944: call PCSetUp(ctx(1),ierr)
1946: call VecDestroy(eold,ierr)
1949: return
1950: end
1951: subroutine Setmat(ctx,ierr)
1953: implicit none
1955: #include petsc/finclude/petsc.h
1957: ! Common blocks:
1958: #include "ex74fcomd.h"
1959: #include "ex74ftube.h"
1961: ! Input/output variables:
1962: PetscFortranAddr ctx(*)
1963: integer ierr
1965: ! Local variables:
1966: PetscScalar lx_v(0:1)
1967: PetscOffset lx_i
1969: double precision xmult, himh, hiph, diag, upper, lower
1970: double precision hi, hip1, him1
1971: double precision rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw
1972: double precision see, sep, sew, seef, sewf, tef, twf, ref, rwf, kef, kwf, xmulte, xmultw
1973: !
1974: integer im, nx, ny
1975: !
1976: ! get pointers to xold
1977: !
1978: call VecGetArray(ctx(4),lx_v,lx_i,ierr)
1979:
1981: !
1982: !############################
1983: !
1984: ! loop over all cells begin
1985: !
1986: !############################
1987: !
1988: do im = 1,mx
1989: !
1990: ! set scalars
1991: !
1992: call Setpbc(im,lx_v(lx_i), rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw)
1993: !
1994: ! set diffusion coefficients
1995: !
1996: see = (erge/rhoe) - (0.5d+0 * ue * ue)
1997: sep = (ergp/rhop) - (0.5d+0 * up * up)
1998: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2000: seef = 0.5d+0 * (see + sep)
2001: sewf = 0.5d+0 * (sew + sep)
2003: tef = seef / csubv
2004: twf = sewf / csubv
2006: ref = 0.5d+0 * (rhoe + rhop)
2007: rwf = 0.5d+0 * (rhow + rhop)
2009: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2010: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2012: if (wilson) then
2013: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2014: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2015: endif
2016: !
2017: ! set coefficients
2018: !
2019: xmult = dt / (dx * dx * csubv)
2021: xmulte = xmult * kef
2022: xmultw = xmult * kwf
2024: upper = -(xmulte / rhoe)
2025: lower = -(xmultw / rhow)
2027: diag = 1.0d+0 + ( (xmulte + xmultw) / rhop )
2029: !
2030: ! load coefficients into the matrix
2031: !
2032: call MatSetValues(ctx(2),1,im-1,1,im-1,diag,INSERT_VALUES,ierr)
2034: if (im .eq. 1) then
2035: call MatSetValues(ctx(2),1,im-1,1,im ,upper, INSERT_VALUES,ierr)
2036: elseif (im .eq. mx) then
2037: call MatSetValues(ctx(2),1,im-1,1,im-2,lower,INSERT_VALUES,ierr)
2038: else
2039: call MatSetValues(ctx(2),1,im-1,1,im ,upper,INSERT_VALUES,ierr)
2040: call MatSetValues(ctx(2),1,im-1,1,im-2,lower, INSERT_VALUES,ierr)
2041: endif
2044: enddo
2045: !
2046: !############################
2047: !
2048: ! loop over all cells end
2049: !
2050: !############################
2051: !
2052:
2053: !
2054: ! final load of matrix
2055: !
2056: call MatAssemblyBegin(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2057: call MatAssemblyEnd(ctx(2),MAT_FINAL_ASSEMBLY,ierr)
2059: if (debug) then
2060: call MatGetSize(ctx(2),nx,ny,ierr)
2061: write(*,*) 'in setup nx = ',nx,' ny = ',ny
2062: call MatView(ctx(2),PETSC_VIEWER_DRAW_WORLD,ierr)
2063: endif
2065: call VecRestoreArray (ctx(4),lx_v,lx_i,ierr)
2069: return
2070: end
2071: subroutine Setpbc(i,x, rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw)
2073: implicit none
2075: ! Common blocks:
2076: #include "ex74fcomd.h"
2078: ! Input/output variables:
2079: PetscScalar x(mx*neq)
2080: integer i
2081: double precision rhoee, rhoe, rhop, rhow, rhoww
2082: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2083: double precision ergee, erge, ergp, ergw, ergww
2084: double precision vele, velp, velw
2086: ! Local variables:
2087: integer jr, jru, je
2089: !
2090: ! set pointers
2091: !
2092: jr = (neq*i) - 2
2093: jru = (neq*i) - 1
2094: je = (neq*i)
2096: if (debug) then
2097: write(*,*)
2098: write(*,*) 'in Setpbc jr,jru,je = ',jr,jru,je
2099: write(*,*)
2100: endif
2101:
2102: if (i .eq. 1) then
2104: rhoee = x(jr+(2*neq))
2105: rhoe = x(jr+neq)
2106: rhop = x(jr)
2107: rhow = x(jr)
2108: rhoww = x(jr)
2110: rhouee = x(jru+(2*neq))
2111: rhoue = x(jru+neq)
2112: rhoup = x(jru)
2113: rhouw = x(jru)
2114: rhouww = x(jru)
2116: ergee = x(je+(2*neq))
2117: erge = x(je+neq)
2118: ergp = x(je)
2119: ergw = x(je)
2120: ergww = x(je)
2122: velw = 0.0d+0
2123: velp = rhoup/rhop
2124: vele = rhoue/rhoe
2126: elseif (i .eq. 2) then
2128: rhoee = x(jr+(2*neq))
2129: rhoe = x(jr+neq)
2130: rhop = x(jr)
2131: rhow = x(jr-neq)
2132: rhoww = x(jr-neq)
2134: rhouee = x(jru+(2*neq))
2135: rhoue = x(jru+neq)
2136: rhoup = x(jru)
2137: rhouw = x(jru-neq)
2138: rhouww = x(jru-neq)
2140: ergee = x(je+(2*neq))
2141: erge = x(je+neq)
2142: ergp = x(je)
2143: ergw = x(je-neq)
2144: ergww = x(je-neq)
2146: velw = rhouw/rhow
2147: velp = rhoup/rhop
2148: vele = rhoue/rhoe
2150: elseif (i .eq. mx-1) then
2152: rhoee = x(jr+neq)
2153: rhoe = x(jr+neq)
2154: rhop = x(jr)
2155: rhow = x(jr-neq)
2156: rhoww = x(jr-(2*neq))
2158: rhouee = x(jru+neq)
2159: rhoue = x(jru+neq)
2160: rhoup = x(jru)
2161: rhouw = x(jru-neq)
2162: rhouww = x(jru-(2*neq))
2164: ergee = x(je+neq)
2165: erge = x(je+neq)
2166: ergp = x(je)
2167: ergw = x(je-neq)
2168: ergww = x(je-(2*neq))
2170: velw = rhouw/rhow
2171: velp = rhoup/rhop
2172: vele = rhoue/rhoe
2174: elseif (i .eq. mx) then
2176: rhoee = x(jr)
2177: rhoe = x(jr)
2178: rhop = x(jr)
2179: rhow = x(jr-neq)
2180: rhoww = x(jr-(2*neq))
2182: rhouee = x(jru)
2183: rhoue = x(jru)
2184: rhoup = x(jru)
2185: rhouw = x(jru-neq)
2186: rhouww = x(jru-(2*neq))
2188: ergee = x(je)
2189: erge = x(je)
2190: ergp = x(je)
2191: ergw = x(je-neq)
2192: ergww = x(je-(2*neq))
2194: velw = rhouw/rhow
2195: velp = rhoup/rhop
2196: vele = 0.0d+0
2198: else
2200: rhoee = x(jr+(2*neq))
2201: rhoe = x(jr+neq)
2202: rhop = x(jr)
2203: rhow = x(jr-neq)
2204: rhoww = x(jr-(2*neq))
2206: rhouee = x(jru+(2*neq))
2207: rhoue = x(jru+neq)
2208: rhoup = x(jru)
2209: rhouw = x(jru-neq)
2210: rhouww = x(jru-(2*neq))
2212: ergee = x(je+(2*neq))
2213: erge = x(je+neq)
2214: ergp = x(je)
2215: ergw = x(je-neq)
2216: ergww = x(je-(2*neq))
2218: velw = rhouw/rhow
2219: velp = rhoup/rhop
2220: vele = rhoue/rhoe
2222: endif
2224: if (debug) then
2225: write(*,*)
2226: write(*,*) 'in Setpbc ',i,jr,jru,je
2227: write(*,*) 'mx = ',mx
2228: write(*,*)
2229: endif
2232: return
2233: end
2234: subroutine Setpbcn(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, ue, up, uw, jbc)
2236: implicit none
2238: ! Common blocks:
2239: #include "ex74fcomd.h"
2241: ! Input/output variables:
2242: integer jbc
2243: double precision rhoee, rhoe, rhop, rhow, rhoww
2244: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2245: double precision ergee, erge, ergp, ergw, ergww
2246: double precision ue, up, uw
2248: ! Local variables:
2250: if (jbc .eq. 1) then
2251: rhoww = rhop
2252: rhow = rhop
2253: rhouww = rhoup
2254: rhouw = rhoup
2255: ergww = ergp
2256: ergw = ergp
2257: uw = 0.0d+0
2258: elseif (jbc .eq. 2) then
2259: rhoww = rhow
2260: rhouww = rhouw
2261: ergww = ergw
2262: uw = rhouw / rhow
2263: else
2264: uw = rhouw / rhow
2265: endif
2267: if (jbc .eq. mx) then
2268: rhoee = rhop
2269: rhoe = rhop
2270: rhouee = rhoup
2271: rhoue = rhoup
2272: ergee = ergp
2273: erge = ergp
2274: ue = 0.0d+0
2275: elseif (jbc .eq. mx-1) then
2276: rhoee = rhoe
2277: rhouee = rhoue
2278: ergee = erge
2279: ue = rhoue / rhoe
2280: else
2281: ue = rhoue / rhoe
2282: endif
2284: up = rhoup / rhop
2286: if (debug) then
2287: write(*,*) 'in Setpbcn ',jbc, 'mx = ',mx
2288: endif
2291: return
2292: end
2293: double precision function cont(rhoee, rhoe, rhop, rhow, rhoww, rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww,jcont,xold)
2294: !
2295: ! This function computes the residual
2296: ! for the 1-D continuity equation
2297: !
2298: !
2299: implicit none
2301: include 'ex74fcomd.h'
2302: include 'ex74ftube.h'
2303: !
2304: ! input variables
2305: !
2306: double precision rhoee, rhoe, rhop, rhow, rhoww
2307: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2308: double precision ergee, erge, ergp, ergw, ergww
2309: double precision xold(mx*neq)
2310: !
2311: integer jcont
2312: !
2313: ! local variables
2314: !
2315: double precision theta1
2316: integer jr
2317: !
2318: ! new
2319: !
2320: double precision velfw, velfe
2321: double precision vele,velp,velw
2322: double precision fluxe, fluxw
2323: double precision urhoe, urhow
2324: double precision source
2325: !
2326: ! old
2327: !
2328: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2329: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2330: double precision ergoee, ergoe, ergop, ergow, ergoww
2331: double precision teoee, teoe, teop, teow, teoww, uoe, uoee, uop, uow, uoww
2332: double precision velfow, velfoe
2333: double precision veloe,velop,velow
2334: double precision fluxoe, fluxow
2335: double precision urhooe, urhoow
2336: double precision sourceo
2337: !
2338: ! functions
2339: !
2340: double precision godunov2
2341: double precision upwind, fluxlim
2342: !
2343: !
2344: ! ******************************************************************
2345: !
2346: !
2347: !
2348: if (debug) then
2349: write(*,*)
2350: write(*,*) 'in cont',jcont,' ihod = ',ihod
2351: write(*,*) 'rhoee = ',rhoee, ' rhoe = ',rhoe
2352: write(*,*) 'rhop = ',rhop
2353: write(*,*) 'rhoww = ',rhoww, ' rhow = ',rhow
2354: write(*,*)
2355: endif
2357: jr = (neq*jcont) - 2
2359: !########################
2360: !
2361: ! NEW
2362: !
2363: !########################
2365: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww, ergee, erge, ergp, ergw, ergww, vele, velp, velw, jcont)
2367: velfe = 0.5d+0 * (vele + velp)
2368: velfw = 0.5d+0 * (velw + velp)
2370: if (ihod .eq. 1) then
2372: urhoe = upwind(rhop,rhoe,velfe)
2373: urhow = upwind(rhow,rhop,velfw)
2375: elseif (ihod .eq. 2) then
2377: urhoe = fluxlim(rhow,rhop,rhoe,rhoee,velfe)
2378: urhow = fluxlim(rhoww,rhow,rhop,rhoe,velfw)
2380: endif
2382: if (ihod .eq. 3) then
2383: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,1)
2384: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,1)
2385: else
2386: fluxe = (dt/dx) * urhoe
2387: fluxw = (dt/dx) * urhow
2388: endif
2391: source = 0.0d+0
2393: !########################
2394: !
2395: ! OLD
2396: !
2397: !########################
2399: call Setpbc(jcont,xold,rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
2401: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jcont)
2403: velfoe = 0.5d+0 * (veloe + velop)
2404: velfow = 0.5d+0 * (velow + velop)
2407: if (ihod .eq. 1) then
2409: urhooe = upwind(rhoop,rhooe,velfoe)
2410: urhoow = upwind(rhoow,rhoop,velfow)
2412: elseif (ihod .eq. 2) then
2414: urhooe = fluxlim(rhoow,rhoop,rhooe,rhooee,velfoe)
2415: urhoow = fluxlim(rhooww,rhoow,rhoop,rhooe,velfow)
2417: endif
2419: if (ihod .eq. 3) then
2420: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,1)
2421: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,1)
2422: else
2423: fluxoe = (dt/dx) * urhooe
2424: fluxow = (dt/dx) * urhoow
2425: endif
2427: sourceo = 0.0d+0
2430: !########################
2431: !
2432: ! FUNCTION
2433: !
2434: !########################
2436: theta1 = 1.0d+0 - theta
2437: cont = (rhop - xold(jr)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2438: !VAM
2439: if (probnum .eq. 3) then
2440: cont = 0.0d+0
2441: endif
2442: !VAM
2445: if (debug) then
2446: write(*,*)
2447: write(*,*) 'cont(',jcont,') = ',cont
2448: write(*,*) 'theta = ',theta,'rhop = ',rhop
2449: write(*,*) 'source = ',source,' sourceo = ',sourceo
2450: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
2451: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
2452: write(*,*) 'urhoe = ',urhoe,' urhow = ',urhow
2453: write(*,*) 'urhooe = ',urhooe,' urhoow = ',urhoow
2454: write(*,*)
2455: endif
2457: return
2458: end
2459: double precision function eexact(x,t)
2461: implicit none
2463: double precision x,t
2464: double precision xot, head, tail, contact, ufan
2465: double precision xpow, grat, urat
2466: double precision uexact
2469: logical debug
2471: include 'ex74ftube.h'
2473: debug = .false.
2476: if (t .le. 0.0d+0) then
2477: if (x .gt. 0.0d+0) then
2478: eexact = e1
2479: else
2480: eexact = e4
2481: endif
2482: else
2484: xot = x/t
2485: head = -a4
2486: tail = v3 - a3
2487: contact = v2
2489: if (xot .lt. head) then
2490: eexact = e4
2491: elseif (xot .gt. sspd) then
2492: eexact = e1
2493: elseif (xot .gt. contact) then
2494: eexact = e2
2495: elseif (xot .gt. tail) then
2496: eexact = e3
2497: else
2498: ufan = uexact(x,t)
2499: grat = (gamma - 1.0d+0) / 2.0d+0
2500: xpow = 2.0d+0
2501: urat = ufan / a4
2502: eexact = e4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
2503: endif
2505: endif
2508: if (debug) then
2509: write(*,*)
2510: write(*,*) 'eexact(',x,',',t,') = ',eexact
2511: write(*,*)
2512: endif
2514: return
2515: end
2516: subroutine eigen(ht,uht)
2517: !23456789012345678901234567890123456789012345678901234567890123456789012
2518: !
2519: ! subroutine eigen
2520: !
2521: ! This subroutine computes the eigen values and eigen vectors
2522: !
2523: !23456789012345678901234567890123456789012345678901234567890123456789012
2526: !#######################################################################
2528: implicit none
2530: include 'ex74fcomd.h'
2532: double precision ht, uht
2534: double precision ut, at, lam1, lam2
2537: !#######################################################################
2539: ut = uht / ht
2540: at = sqrt( ht)
2542: lam1 = ut - at
2543: lam2 = ut + at
2545: eigval(1) = lam1
2546: eigval(2) = lam2
2548: eigvec(1,1) = 1.0d+0
2549: eigvec(2,1) = lam1
2550: eigvec(1,2) = 1.0d+0
2551: eigvec(2,2) = lam2
2553: rinv(1,1) = lam2 / (2.0d+0 * at)
2554: rinv(2,1) = -lam1 / (2.0d+0 * at)
2555: rinv(1,2) = -1.0d+0 / (2.0d+0 * at)
2556: rinv(2,2) = 1.0d+0 / (2.0d+0 * at)
2559: return
2560: end
2561: subroutine eigene(r,ru,e,l1,l2,l3)
2562: !23456789012345678901234567890123456789012345678901234567890123456789012
2563: !
2564: ! subroutine eigene
2565: !
2566: ! This subroutine computes the eigen values for the entropy fix
2567: !
2568: !23456789012345678901234567890123456789012345678901234567890123456789012
2571: !#######################################################################
2573: implicit none
2575: include 'ex74ftube.h'
2577: double precision r,ru,e,l1,l2,l3
2579: double precision p,u,a
2581: double precision eos
2582: integer ierr
2584: logical debug
2587: !#######################################################################
2589: debug = .false.
2591: if (debug) then
2592: write(*,*)
2593: write(*,*) 'gamma = ',gamma
2594: write(*,*) 'r,ru,e = ',r,ru,e
2595: write(*,*)
2596: endif
2598: p = eos(r,ru,e)
2599: u = ru/r
2600: if ( ((gamma * p)/r) .lt. 0.0d+0 ) then
2601: write(*,*)
2602: write(*,*) 'gamma = ',gamma
2603: write(*,*) 'r = ',r
2604: write(*,*) 'p = ',p
2605: write(*,*)
2606: call PetscFinalize(ierr)
2607: stop
2608: endif
2609: a = sqrt((gamma * p)/r)
2611: if (debug) then
2612: write(*,*)
2613: write(*,*) 'p,u,a = ',p,u,a
2614: write(*,*)
2615: endif
2617: l1 = u - a
2618: l2 = u
2619: l3 = u + a
2621: return
2622: end
2623: double precision function energy(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,jerg,xold)
2624: !
2625: ! This function computes the residual
2626: ! for the 1-D energy equation
2627: !
2628: !
2629: implicit none
2631: include 'ex74fcomd.h'
2632: include 'ex74ftube.h'
2633: !
2634: ! input variables
2635: !
2636: double precision rhoee, rhoe, rhop, rhow, rhoww
2637: double precision rhouee, rhoue, rhoup, rhouw, rhouww
2638: double precision ergee, erge, ergp, ergw, ergww
2639: double precision xold(mx*neq)
2640: !
2641: integer jerg
2642: !
2643: ! local variables
2644: !
2645: double precision theta1
2646: integer je
2647: !
2648: ! new
2649: !
2650: double precision velfw, velfe
2651: double precision vele,velp,velw
2652: double precision fluxe, fluxw
2653: double precision uepe, uepw
2654: double precision ue, up, uw
2655: double precision see, sep, sew
2656: double precision seef, sewf
2657: double precision upe, upw
2658: double precision presse, pressw
2659: double precision source
2660: double precision te, tp, tw
2661: double precision tef, twf, ref, rwf
2662: double precision kef, kwf
2663: double precision hflxe, hflxw
2664: !
2665: ! old
2666: !
2667: double precision rhooee, rhooe, rhoop, rhoow, rhooww
2668: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
2669: double precision ergoee, ergoe, ergop, ergow, ergoww
2670: double precision velfow, velfoe
2671: double precision veloe,velop,velow
2672: double precision fluxoe, fluxow
2673: double precision uepoe, uepow
2674: double precision uoe, uop, uow
2675: double precision seoe, seop, seow
2676: double precision seoef, seowf
2677: double precision upoe, upow
2678: double precision pressoe, pressow
2679: double precision sourceo
2680: double precision toe, top, tow
2681: double precision toef, towf, roef, rowf
2682: double precision koef, kowf
2683: double precision hflxoe, hflxow
2684: !
2685: ! functions
2686: !
2687: double precision godunov2, eos
2688: double precision upwind, fluxlim
2690: !
2691: !
2692: ! ******************************************************************
2693: !
2694: !
2695: !
2696: je = (neq*jerg)
2698: !########################
2699: !
2700: ! NEW
2701: !
2702: !########################
2704: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,vele, velp, velw, jerg)
2706: pressw = eos(rhow, rhouw, ergw)
2707: presse = eos(rhoe, rhoue, erge)
2709: uw = rhouw / rhow
2710: up = rhoup / rhop
2711: ue = rhoue / rhoe
2713: upw = uw * pressw
2714: upe = ue * presse
2716: velfe = 0.5d+0 * (vele + velp)
2717: velfw = 0.5d+0 * (velw + velp)
2719: if (ihod .eq. 1) then
2721: uepe = upwind(ergp,erge,velfe)
2722: uepw = upwind(ergw,ergp,velfw)
2724: elseif (ihod .eq. 2) then
2726: uepe = fluxlim(ergw,ergp,erge,ergee,velfe)
2727: uepw = fluxlim(ergww,ergw,ergp,erge,velfw)
2729: endif
2731: if (ihod .eq. 3) then
2732: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee, rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,3)
2733: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,3)
2734: else
2735: fluxe = (dt/dx) * ( uepe + (0.5d+0*upe) )
2736: fluxw = (dt/dx) * ( uepw + (0.5d+0*upw) )
2737: endif
2738: !
2739: ! radiation
2740: !
2741: if (kappa0 .eq. 0.0d+0) then
2742: source = 0.0d+0
2743: else
2745: see = (erge/rhoe) - (0.5d+0 * ue * ue)
2746: sep = (ergp/rhop) - (0.5d+0 * up * up)
2747: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
2749: seef = 0.5d+0 * (see + sep)
2750: sewf = 0.5d+0 * (sew + sep)
2752: te = see / csubv
2753: tp = sep / csubv
2754: tw = sew / csubv
2756: tef = seef / csubv
2757: twf = sewf / csubv
2759: ref = 0.5d+0 * (rhoe + rhop)
2760: rwf = 0.5d+0 * (rhow + rhop)
2762: kef = kappa0 * (ref ** kappaa) * (tef ** kappab)
2763: kwf = kappa0 * (rwf ** kappaa) * (twf ** kappab)
2765: if (wilson) then
2766: kef = 1.0d+0 / ((1.0d+0/kef)+(abs(see-sep)/(seef*dx)))
2767: kwf = 1.0d+0 / ((1.0d+0/kwf)+(abs(sep-sew)/(sewf*dx)))
2768: endif
2770: if ( debug .and. (kef .gt. 1.0d+10) ) then
2771: write(*,*) 'kef = ',kef,ref,tef,kappaa,kappab,kappa0
2772: endif
2773: if ( debug .and. (kwf .gt. 1.0d+10) ) then
2774: write(*,*) 'kwf = ',kwf,rwf,twf,kappaa,kappab,kappa0
2775: endif
2777: hflxe = kef * (te - tp) / dx
2778: hflxw = kwf * (tp - tw) / dx
2780: source = (dt/dx) * (hflxe - hflxw)
2782: endif
2784: !########################
2785: !
2786: ! OLD
2787: !
2788: !########################
2790: call Setpbc(jerg,xold, rhooee, rhooe, rhoop, rhoow, rhooww,rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
2792: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jerg)
2794: pressow = eos(rhoow, rhouow, ergow)
2795: pressoe = eos(rhooe, rhouoe, ergoe)
2797: uow = rhouow / rhoow
2798: uop = rhouop / rhoop
2799: uoe = rhouoe / rhooe
2801: upow = uow * pressow
2802: upoe = uoe * pressoe
2804: velfoe = 0.5d+0 * (veloe + velop)
2805: velfow = 0.5d+0 * (velow + velop)
2808: if (ihod .eq. 1) then
2810: uepoe = upwind(ergop,ergoe,velfoe)
2811: uepow = upwind(ergow,ergop,velfow)
2813: elseif (ihod .eq. 2) then
2815: uepoe = fluxlim(ergow,ergop,ergoe,ergoee,velfoe)
2816: uepow = fluxlim(ergoww,ergow,ergop,ergoe,velfow)
2818: endif
2820: if (ihod .eq. 3) then
2821: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee,rhouow,rhouop,rhouoe,rhouoee, ergow, ergop, ergoe, ergoee,3)
2822: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe, rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,3)
2823: else
2824: fluxoe = (dt/dx) * ( uepoe + (0.5d+0 * upoe) )
2825: fluxow = (dt/dx) * ( uepow + (0.5d+0 * upow) )
2826: endif
2828: !
2829: ! old radiation
2830: !
2831: if (kappa0 .eq. 0.0d+0) then
2832: sourceo = 0.0d+0
2833: else
2835: seoe = (ergoe/rhooe) - (0.5d+0 * uoe * uoe)
2836: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
2837: seow = (ergow/rhoow) - (0.5d+0 * uow * uow)
2839: seoef = 0.5d+0 * (seoe + seop)
2840: seowf = 0.5d+0 * (seow + seop)
2842: toe = seoe / csubv
2843: top = seop / csubv
2844: tow = seow / csubv
2846: toef = seoef / csubv
2847: towf = seowf / csubv
2849: roef = 0.5d+0 * (rhooe + rhoop)
2850: rowf = 0.5d+0 * (rhoow + rhoop)
2852: koef = kappa0 * (roef ** kappaa) * (toef ** kappab)
2853: kowf = kappa0 * (rowf ** kappaa) * (towf ** kappab)
2855: if (wilson) then
2856: koef = 1.0d+0 / ((1.0d+0/koef)+(abs(seoe-seop)/(seoef*dx)))
2857: kowf = 1.0d+0 / ((1.0d+0/kowf)+(abs(seop-seow)/(seowf*dx)))
2858: endif
2860: if ( debug .and. (koef .gt. 1.0d+10) ) then
2861: write(*,*) 'koef = ',koef,roef,toef,kappaa,kappab,kappa0
2862: endif
2863: if ( debug .and. (kowf .gt. 1.0d+10) ) then
2864: write(*,*) 'kowf = ',kowf,rowf,towf,kappaa,kappab,kappa0
2865: endif
2867: hflxoe = koef * (toe - top) / dx
2868: hflxow = kowf * (top - tow) / dx
2870: sourceo = (dt/dx) * (hflxoe - hflxow)
2872: endif
2875: !########################
2876: !
2877: ! FUNCTION
2878: !
2879: !########################
2881: !VAM
2882: if (probnum .eq. 3) then
2883: fluxe = 0.0d+0
2884: fluxw = 0.0d+0
2885: fluxoe = 0.0d+0
2886: fluxow = 0.0d+0
2887: endif
2888: !VAM
2890: theta1 = 1.0d+0 - theta
2891: energy = (ergp - xold(je)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
2893: if (debug) then
2894: write(*,*)
2895: write(*,*) 'energy(',jerg,') = ',energy
2896: write(*,*)
2897: write(*,*) fluxe,fluxw
2898: write(*,*) fluxoe,fluxow
2899: write(*,*) source,sourceo
2900: write(*,*)
2901: endif
2903: return
2904: end
2905: double precision function eos(r,ru,e)
2907: implicit none
2909: double precision r,ru,e
2911: double precision se, u
2913: integer ierr
2915: logical debug
2917: include "ex74ftube.h"
2919: debug = .false.
2921: if (debug) then
2922: write(*,*)
2923: write(*,*) 'in eos r,ru,e'
2924: write(*,*) r,ru,e
2925: write(*,*)
2926: endif
2928: u = ru/r
2930: se = (e/r) - (0.5d+0 * u * u)
2931: eos = (gamma - 1.0d+0) * r * se
2933: if (eos .lt. 0.0d+0) then
2934: write(*,*)
2935: write(*,*) 'eos = ',eos
2936: write(*,*) 'gamma = ',gamma
2937: write(*,*) 'r = ',r
2938: write(*,*) 'se = ',se
2939: write(*,*) 'e = ',e
2940: write(*,*) 'u = ',u
2941: write(*,*) 'ru = ',ru
2942: call PetscFinalize(ierr)
2943: write(*,*)
2944: stop
2945: endif
2947: if (debug) then
2948: write(*,*)
2949: write(*,*) 'in eos u,se,eos'
2950: write(*,*) u,se,eos
2951: write(*,*)
2952: endif
2955: return
2956: end
2957: subroutine eval2
2959: implicit none
2961: double precision prat, grat, xnum, xdenom
2964: logical debug
2966: include 'ex74ftube.h'
2968: debug = .false.
2970: prat = p2/p1
2971: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
2973: xnum = grat + prat
2974: xdenom = 1.0d+0 + (prat * grat)
2975:
2976: e2 = e1 * prat * (xnum/xdenom)
2977:
2980: if (debug) then
2981: write(*,*)
2982: write(*,*) 'e1 = ',e1
2983: write(*,*) 'e2 = ',e2
2984: endif
2986: return
2987: end
2988: subroutine exact0
2990: implicit none
2992: double precision tol, xn
2993: double precision shockp, fprime
2995: integer maxnewt, niter
2997: logical found, debug
2999: include 'ex74ftube.h'
3001: debug = .false.
3003: tol = 1.0d-10
3005: maxnewt = 40
3006:
3007: a1 = sqrt(gamma*p1/r1)
3008: a4 = sqrt(gamma*p4/r4)
3012: found = .false.
3013: niter = 0
3015: xn = 0.5d+0 * (p1 + p4)
3017: 10 if ( (.not. found) .and. (niter .le. maxnewt) ) then
3019: niter = niter + 1
3021: xn = xn - (shockp(xn) / fprime(xn))
3023: if (debug) then
3024: write(*,*) niter,shockp(xn),xn
3025: endif
3027: if ( abs(shockp(xn)) .lt. tol ) then
3028: found = .true.
3029: endif
3031: goto 10
3033: endif
3035: if (.not. found) then
3037: write(*,*) 'newton failed'
3038: write(*,*) xn,shockp(xn)
3039: stop
3041: endif
3043: p2 = xn
3046: if (debug) then
3047: write(*,*)
3048: write(*,*) 'p1 = ',p1
3049: write(*,*) 'p2 = ',p2
3050: write(*,*) 'p4 = ',p4
3051: write(*,*)
3052: endif
3054: return
3055: end
3056: double precision function flux(r,ru,e,eqn)
3057: !23456789012345678901234567890123456789012345678901234567890123456789012
3058: !
3059: ! function flux
3060: !
3061: ! This function computes the flux at a face
3062: !
3063: !23456789012345678901234567890123456789012345678901234567890123456789012
3066: !#######################################################################
3068: implicit none
3070: include 'ex74fcomd.h'
3071: include 'ex74ftube.h'
3073: double precision r, ru, e
3075: integer eqn
3077: double precision p,u
3079: double precision eos
3082: !#######################################################################
3084: p = eos(r,ru,e)
3085: u = ru/r
3087: if (eqn .eq. 1) then
3088: flux = ru
3089: elseif (eqn .eq. 2) then
3090: flux = (u * ru) + p
3091: else
3092: flux = u * (e + p)
3093: endif
3095: return
3096: end
3097: double precision function fluxlim(fww,fw,fe,fee,vp)
3098: !23456789012345678901234567890123456789012345678901234567890123456789012
3099: !
3100: ! function fluxlim
3101: !
3102: ! this function computes the flux limited quick face value
3103: !
3104: !23456789012345678901234567890123456789012345678901234567890123456789012
3107: !#######################################################################
3109: implicit none
3111: double precision fww, fw, fe, fee, vp
3113: double precision fd, fc, fu
3115: double precision f1, f2, f3, f4, fhod, beta, flc
3117: double precision med, quick
3118:
3119: logical limit
3121: !#######################################################################
3123: limit = .true.
3125: if (vp .gt. 0.0d+0) then
3126: fd = fe
3127: fc = fw
3128: fu = fww
3129: else
3130: fd = fw
3131: fc = fe
3132: fu = fee
3133: endif
3135: fhod = quick(fd,fc,fu)
3137: if (limit) then
3139: beta = 0.25d+0
3140: flc = 4.0d+0
3142: f1 = fc
3143: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3144: f3 = fu + ( flc * (fc - fu) )
3145: f4 = med(f1,f2,f3)
3146: fluxlim = vp * med(f1,f4,fhod)
3148: else
3150: fluxlim = vp * fhod
3152: endif
3154: return
3155: end
3156: double precision function fluxlim2(fww,fw,fe,fee,vp)
3157: !23456789012345678901234567890123456789012345678901234567890123456789012
3158: !
3159: ! function fluxlim2
3160: !
3161: ! this function computes the flux limited quick face value
3162: !
3163: !23456789012345678901234567890123456789012345678901234567890123456789012
3166: !#######################################################################
3168: implicit none
3170: double precision fww, fw, fe, fee, vp
3172: double precision fd, fc, fu
3174: double precision f1, f2, f3, f4, fhod, beta, flc
3176: double precision med, quick
3177:
3178: logical limit, debug
3180: !#######################################################################
3182: debug = .false.
3184: if (debug) then
3185: write(*,*)
3186: write(*,*) 'in fluxlim2 fee,fe,fw,fww'
3187: write(*,*) fee,fe,fw,fww
3188: write(*,*)
3189: endif
3191: limit = .true.
3193: if (vp .gt. 0.0d+0) then
3194: fd = fe
3195: fc = fw
3196: fu = fww
3197: else
3198: fd = fw
3199: fc = fe
3200: fu = fee
3201: endif
3203: fhod = quick(fd,fc,fu)
3205: if (limit) then
3207: beta = 0.25d+0
3208: flc = 4.0d+0
3210: f1 = fc
3211: f2 = (beta*fc) + ( (1.0d+0-beta)*fd )
3212: f3 = fu + ( flc * (fc - fu) )
3213: f4 = med(f1,f2,f3)
3214: fluxlim2 = med(f1,f4,fhod)
3216: else
3218: fluxlim2 = fhod
3220: endif
3222: return
3223: end
3224: double precision function fprime(x)
3226: implicit none
3228: double precision x, eps
3229: double precision shockp
3231: eps = 1.0d-8
3233: fprime = ( shockp(x+eps) - shockp(x) ) / eps
3235: return
3236: end
3237: double precision function godunov2(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,eqn)
3238: !23456789012345678901234567890123456789012345678901234567890123456789012
3239: !
3240: ! function godunov2
3241: !
3242: ! this function computes the roe/godunov2 face value
3243: !
3244: !23456789012345678901234567890123456789012345678901234567890123456789012
3247: !#######################################################################
3249: implicit none
3251: include 'ex74fcomd.h'
3252: include 'ex74ftube.h'
3254: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3256: integer eqn
3258: double precision rrg, rlg, rurg, rulg, erg, elg
3260: double precision hlle
3263: !#######################################################################
3265: if (gorder .eq. 1) then
3266: rrg = rr
3267: rlg = rl
3268: rurg = rur
3269: rulg = rul
3270: erg = er
3271: elg = el
3272: else
3273: call secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err,rrg, rlg,rurg, rulg, erg, elg)
3274: endif
3276: !VAM if (ientro .eq. 0) then
3277: !VAM godunov2 = godunov(uhlg,uhrg,hlg,hrg,eqn)
3278: !VAM elseif(ientro .eq. 1) then
3279: !VAM godunov2 = godent(uhlg,uhrg,hlg,hrg,eqn)
3280: !VAM else
3281: godunov2 = hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3282: !VAM endif
3285: return
3286: end
3287: double precision function hlle(rrg,rlg,rurg,rulg,erg,elg,eqn)
3288: !23456789012345678901234567890123456789012345678901234567890123456789012
3289: !
3290: ! function hlle
3291: !
3292: ! this function computes the roe/hlle face value
3293: !
3294: !23456789012345678901234567890123456789012345678901234567890123456789012
3297: !#######################################################################
3299: implicit none
3301: include 'ex74fcomd.h'
3302: include 'ex74ftube.h'
3304: double precision rrg,rlg,rurg,rulg,erg,elg
3305: integer eqn
3307: double precision laml1, laml2, laml3
3308: double precision lamr1, lamr2, lamr3
3309: double precision sl, sr
3312: double precision flux
3314: integer i, j, ispeed
3315:
3317: !#######################################################################
3319: ispeed = 1
3321: do i = 1,neq
3322: fr(i) = flux(rrg,rurg,erg,i)
3323: fl(i) = flux(rlg,rulg,elg,i)
3324: enddo
3326: deltau(1) = rrg - rlg
3327: deltau(2) = rurg - rulg
3328: deltau(3) = erg - elg
3330: !VAM call roestat(uhl,uhr,hl,hr,ht,uht)
3332: !VAM call eigene(ht,uht,lamt1, lamt2)
3333: call eigene(rrg,rurg,erg,lamr1,lamr2,lamr3)
3334: call eigene(rlg,rulg,elg,laml1,laml2,laml3)
3336: !VAM if (ispeed .eq. 1) then
3337: !VAM sl = min(laml1,lamt1)
3338: !VAM sr = max(lamt2,lamr2)
3339: !VAM else
3340: sl = min(laml1,lamr1)
3341: sr = max(laml3,lamr3)
3342: !VAM endif
3345: do i = 1,neq
3346: froe(i) = ( (sr*fl(i)) - (sl*fr(i)) + (sl*sr*deltau(i)) )/(sr-sl)
3347: enddo
3349: hlle = froe(eqn)
3352: return
3353: end
3354: double precision function med(x1,x2,x3)
3355: !23456789012345678901234567890123456789012345678901234567890123456789012
3356: !
3357: ! function med
3358: !
3359: ! this function computes the median of three numbers
3360: !
3361: !23456789012345678901234567890123456789012345678901234567890123456789012
3364: !#######################################################################
3366: implicit none
3368: double precision x1, x2, x3
3369: double precision xhi, xlo
3371: !#######################################################################
3373: xhi = max(x1,x2,x3)
3374: xlo = min(x1,x2,x3)
3376: med = x1 + x2 + x3 - xhi - xlo
3378: return
3379: end
3380: double precision function mom(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,jmom,xold)
3381: !
3382: ! This function computes the residual
3383: ! for the 1-D momentum equation
3384: !
3385: !
3386: implicit none
3388: include 'ex74fcomd.h'
3389: include 'ex74ftube.h'
3390: !
3391: ! input variables
3392: !
3393: double precision rhoee, rhoe, rhop, rhow, rhoww
3394: double precision rhouee, rhoue, rhoup, rhouw, rhouww
3395: double precision ergee, erge, ergp, ergw, ergww
3396: double precision xold(mx*neq)
3397: !
3398: integer jmom
3399: !
3400: ! local variables
3401: !
3402: double precision theta1
3403: integer jru
3404: !
3405: ! new
3406: !
3407: double precision velfw, velfe
3408: double precision vele,velp,velw
3409: double precision fluxe, fluxw
3410: double precision uurhoe, uurhow
3411: double precision pressee, presse, pressp,pressw, pressww
3412: double precision rupee, rupe, rupp, rupw, rupww
3413: double precision uee, ue, up, uw, uww
3414: double precision source
3415: !
3416: ! old
3417: !
3418: double precision rhooee, rhooe, rhoop, rhoow, rhooww
3419: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
3420: double precision ergoee, ergoe, ergop, ergow, ergoww
3421: double precision velfow, velfoe
3422: double precision veloe,velop,velow
3423: double precision fluxoe, fluxow
3424: double precision uurhooe, uurhoow
3425: double precision pressoee, pressoe, pressop, pressow, pressoww
3426: double precision rupoee, rupoe, rupop, rupow, rupoww
3427: double precision uoee, uoe, uop, uow, uoww
3428: double precision sourceo
3430: double precision eps
3431: !
3432: ! functions
3433: !
3434: double precision godunov2, eos
3435: double precision upwind, fluxlim
3436: !
3437: !
3438: ! ******************************************************************
3439: !
3440: !
3441: eps = 1.0d-32
3442: !
3443: jru = (neq*jmom) - 1
3445: !########################
3446: !
3447: ! NEW
3448: !
3449: !########################
3451: call Setpbcn(rhoee, rhoe, rhop, rhow, rhoww,rhouee, rhoue, rhoup, rhouw, rhouww,ergee, erge, ergp, ergw, ergww,vele, velp, velw, jmom)
3453: presse = eos(rhoe, rhoue, erge )
3454: pressw = eos(rhow, rhouw, ergw )
3456: velfe = 0.5d+0 * (vele + velp)
3457: velfw = 0.5d+0 * (velw + velp)
3459: if (ihod .eq. 1) then
3461: uurhoe = upwind(rhoup,rhoue,velfe)
3462: uurhow = upwind(rhouw,rhoup,velfw)
3464: elseif (ihod .eq. 2) then
3466: uurhoe = fluxlim(rhouw,rhoup,rhoue,rhouee,velfe)
3467: uurhow = fluxlim(rhouww,rhouw,rhoup,rhoue,velfw)
3469: endif
3471: if (ihod .eq. 3) then
3472: fluxe = (dt/dx) * godunov2(rhow, rhop, rhoe, rhoee,rhouw,rhoup,rhoue,rhouee, ergw, ergp, erge, ergee,2)
3473: fluxw = (dt/dx) * godunov2(rhoww, rhow, rhop, rhoe, rhouww,rhouw,rhoup,rhoue, ergww, ergw, ergp, erge,2)
3474: else
3475: fluxe = (dt/dx) * ( uurhoe + (0.5d+0 * presse) )
3476: fluxw = (dt/dx) * ( uurhow + (0.5d+0 * pressw) )
3477: endif
3480: source = 0.0d+0
3482: !########################
3483: !
3484: ! OLD
3485: !
3486: !########################
3488: call Setpbc(jmom,xold,rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww,ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow)
3490: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww, rhouoee, rhouoe, rhouop, rhouow, rhouoww, ergoee, ergoe, ergop, ergow, ergoww, veloe, velop, velow, jmom)
3492: pressoe = eos(rhooe, rhouoe, ergoe)
3493: pressow = eos(rhoow, rhouow, ergow)
3495: velfoe = 0.5d+0 * (veloe + velop)
3496: velfow = 0.5d+0 * (velow + velop)
3498: if (ihod .eq. 1) then
3500: uurhooe = upwind(rhouop,rhouoe,velfoe)
3501: uurhoow = upwind(rhouow,rhouop,velfow)
3503: elseif (ihod .eq. 2) then
3505: uurhooe = fluxlim(rhouow,rhouop,rhouoe,rhouoee,velfoe)
3506: uurhoow = fluxlim(rhouoww,rhouow,rhouop,rhouoe,velfow)
3508: endif
3510: if (ihod .eq. 3) then
3511: fluxoe = (dt/dx) * godunov2(rhoow, rhoop, rhooe, rhooee, rhouow,rhouop,rhouoe,rhouoee,ergow, ergop, ergoe, ergoee,2)
3512: fluxow = (dt/dx) * godunov2(rhooww, rhoow, rhoop, rhooe,rhouoww,rhouow,rhouop,rhouoe, ergoww, ergow, ergop, ergoe,2)
3513: else
3514: fluxoe = (dt/dx) * ( uurhooe + (0.5d+0 * pressoe) )
3515: fluxow = (dt/dx) * ( uurhoow + (0.5d+0 * pressow) )
3516: endif
3518: sourceo = 0.0d+0
3521: !########################
3522: !
3523: ! FUNCTION
3524: !
3525: !########################
3527: theta1 = 1.0d+0 - theta
3528: mom = (rhoup - xold(jru)) + ( theta * ( (fluxe - fluxw ) - source ) ) + ( theta1 * ( (fluxoe - fluxow) - sourceo ) )
3529: !VAM
3530: if (probnum .eq. 3) then
3531: mom = 0.0d+0
3532: endif
3533: !VAM
3534: if (debug) then
3535: write(*,*)
3536: write(*,*) 'mom(',jmom,') = ',mom,' theta = ',theta
3537: write(*,*) 'fluxe = ',fluxe,' fluxw = ',fluxw
3538: write(*,*) 'fluxoe = ',fluxoe,' fluxow = ',fluxow
3539: write(*,*) 'presse = ',presse,'pressw = ',pressw
3540: write(*,*) 'pressoe = ',pressoe,'pressow = ',pressow
3541: write(*,*)
3542: endif
3544: return
3545: end
3546: double precision function quick(fd, fc, fu)
3547: !23456789012345678901234567890123456789012345678901234567890123456789012
3548: !
3549: ! function quick
3550: !
3551: ! this function computes the quick face value
3552: !
3553: !23456789012345678901234567890123456789012345678901234567890123456789012
3556: !#######################################################################
3558: implicit none
3560: double precision fd, fc, fu
3562: !#######################################################################
3564: quick = ( (3.0d+0 * fd) + (6.0d+0 * fc) - fu ) / 8.0d+0
3566: return
3567: end
3568: double precision function rexact(x,t)
3570: implicit none
3572: double precision x,t
3573: double precision xot, head, tail, contact, ufan
3574: double precision xpow, grat, urat
3575: double precision uexact
3578: logical debug
3580: include 'ex74ftube.h'
3582: debug = .false.
3585: if (t .le. 0.0d+0) then
3586: if (x .gt. 0.0d+0) then
3587: rexact = r1
3588: else
3589: rexact = r4
3590: endif
3591: else
3593: xot = x/t
3594: head = -a4
3595: tail = v3 - a3
3596: contact = v2
3598: if (xot .lt. head) then
3599: rexact = r4
3600: elseif (xot .gt. sspd) then
3601: rexact = r1
3602: elseif (xot .gt. contact) then
3603: rexact = r2
3604: elseif (xot .gt. tail) then
3605: rexact = r3
3606: else
3607: ufan = uexact(x,t)
3608: grat = (gamma - 1.0d+0) / 2.0d+0
3609: xpow = 1.0d+0 / grat
3610: urat = ufan / a4
3611: rexact = r4 * ( ( 1.0d+0 - (grat * urat) ) ** xpow )
3612: endif
3614: endif
3617: if (debug) then
3618: write(*,*)
3619: write(*,*) 'rexact(',x,',',t,') = ',rexact
3620: write(*,*)
3621: endif
3623: return
3624: end
3625: subroutine roestat(uhl, uhr, hl,hr,ht,uht)
3626: !23456789012345678901234567890123456789012345678901234567890123456789012
3627: !
3628: ! subroutine roestat
3629: !
3630: ! This subroutine computes the roe state at a face
3631: !
3632: !23456789012345678901234567890123456789012345678901234567890123456789012
3635: !#######################################################################
3637: implicit none
3639: include 'ex74fcomd.h'
3641: double precision uhl, uhr, hl, hr, ht, uht
3643: double precision ul, ur, shl, shr, xnum, xdenom
3644:
3647: !#######################################################################
3649: ul = uhl / hl
3650: ur = uhr / hr
3652: shl = sqrt(hl)
3653: shr = sqrt(hr)
3655: xnum = (shl * ul) + (shr * ur)
3656: xdenom = shl + shr
3658: ht = 0.5d+0 * (hl + hr)
3659: uht = ht * ( xnum / xdenom )
3661: return
3662: end
3663: subroutine rval2
3665: implicit none
3667: double precision prat, grat, xnum, xdenom
3670: logical debug
3672: include 'ex74ftube.h'
3674: debug = .false.
3676: prat = p2/p1
3677: grat = (gamma + 1.0d+0) / (gamma - 1.0d+0)
3679: xnum = 1.0d+0 + (grat * prat)
3680: xdenom = grat + prat
3681:
3682: r2 = r1 * (xnum/xdenom)
3683:
3686: if (debug) then
3687: write(*,*)
3688: write(*,*) 'r1 = ',r1
3689: write(*,*) 'r2 = ',r2
3690: endif
3692: return
3693: end
3694: subroutine secondq(rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err, rrg, rlg,rurg, rulg, erg, elg)
3695: !23456789012345678901234567890123456789012345678901234567890123456789012
3696: !
3697: ! subroutine secondq
3698: !
3699: ! this subroutine computes the second order (based on quick) left
3700: ! and right states for the godunov solver.
3701: !
3702: !23456789012345678901234567890123456789012345678901234567890123456789012
3705: !#######################################################################
3707: implicit none
3709: include 'ex74fcomd.h'
3711: double precision rll,rl,rr,rrr,rull,rul,rur,rurr,ell,el,er,err
3712: double precision rrg, rlg,rurg, rulg, erg, elg
3716: double precision veld, ull,ul,ur,urr, ulg, urg
3718: double precision fluxlim2
3721: !#######################################################################
3723: !
3724: ! compute the velocities
3725: !
3726: ull = rull/rll
3727: ul = rul /rl
3728: ur = rur /rr
3729: urr = rurr/rrr
3731: !
3732: ! compute the left state first
3733: !
3734: veld = 1.0d+0
3736: rlg = fluxlim2(rll,rl,rr,rrr,veld)
3737: ulg = fluxlim2(ull,ul,ur,urr,veld)
3738: rulg = rlg * ulg
3739: elg = fluxlim2(ell,el,er,err,veld)
3740: !
3741: ! now compute the right state
3742: !
3743: veld = -1.0d+0
3745: rrg = fluxlim2(rll,rl,rr,rrr,veld)
3746: urg = fluxlim2(ull,ul,ur,urr,veld)
3747: rurg = rrg * urg
3748: erg = fluxlim2(ell,el,er,err,veld)
3752: return
3753: end
3754: double precision function shockp(x)
3756: implicit none
3758: double precision x
3759: double precision xnum, xdenom, xpow, prat, prat2, prat4, gm, gp
3760: logical debug
3762: include 'ex74ftube.h'
3764: debug = .false.
3767: if (debug) then
3768: write(*,*)
3769: write(*,*) 'gamma = ',gamma
3770: write(*,*) 'a1 = ',a1
3771: write(*,*) 'a4 = ',a4
3772: write(*,*) 'p1 = ',p1
3773: write(*,*) 'p2 = ',x
3774: write(*,*)
3775: endif
3777: xnum = (gamma - 1.0d+0) * (a1/a4) * ( (x/p1) - 1.0d+0 )
3778: xdenom = sqrt ( 2.0d+0 * gamma * ( (2.0d+0*gamma) + (gamma + 1.0d+0) * ((x/p1) - 1) ) )
3779: xpow = (-2.0d+0 * gamma) / (gamma - 1.0d+0)
3781: shockp = (x/p1)*((1.0d+0-(xnum/xdenom))**xpow) - (p4/p1)
3784: if (debug) then
3785: write(*,*)
3786: write(*,*) 'xnum = ',xnum
3787: write(*,*) 'gamma = ',gamma
3788: write(*,*) 'a1 = ',a1
3789: write(*,*) 'a4 = ',a4
3790: write(*,*) 'p1 = ',p1
3791: write(*,*) 'xdenom = ',xdenom
3792: write(*,*) 'xpow = ',xpow
3793: write(*,*) 'shockp = ',shockp
3794: write(*,*) 'p2 = ',x
3795: write(*,*)
3796: endif
3798: return
3799: end
3800: double precision function uexact(x,t)
3802: implicit none
3804: double precision x,t
3805: double precision xot, head, tail
3808: logical debug
3810: include 'ex74ftube.h'
3812: debug = .false.
3814: if (debug) then
3815: write(*,*)
3816: write(*,*) 't = ',t
3817: write(*,*) 'x = ',x
3818: write(*,*) 'a4 = ',a4
3819: write(*,*) 'v3 = ',v3
3820: write(*,*) 'a3 = ',a3
3821: write(*,*)
3822: endif
3824: if (t .le. 0.0d+0) then
3825: uexact = 0.0d+0
3826: else
3828: xot = x/t
3829: head = -a4
3830: tail = v3 - a3
3832: if (xot .lt. head) then
3833: uexact = 0.0d+0
3834: elseif (xot .gt. sspd) then
3835: uexact = 0.0d+0
3836: elseif (xot .gt. tail) then
3837: uexact = v2
3838: else
3839: uexact = (2.0d+0 / (gamma + 1.0d+0))* (a4 + xot)
3840: endif
3842: endif
3845: if (debug) then
3846: write(*,*)
3847: !VAM write(*,*) 'x = ',x,' t = ',t
3848: write(*,*) 'uexact = ',uexact
3849: write(*,*)
3850: endif
3852: return
3853: end
3854: double precision function upwind(fw, fe, vp)
3855: !23456789012345678901234567890123456789012345678901234567890123456789012
3856: !
3857: ! function upwind
3858: !
3859: ! this function computes the upwind face value
3860: !
3861: !23456789012345678901234567890123456789012345678901234567890123456789012
3864: !#######################################################################
3866: implicit none
3868: double precision fw, fe, vp
3870: !#######################################################################
3872: if (vp .gt. 0.0) then
3873: upwind = vp * fw
3874: else
3875: upwind = vp * fe
3876: endif
3878: return
3879: end
3880: subroutine uval2
3882: implicit none
3884: double precision prat, grat1, grat2, arat, xnum
3887: logical debug
3889: include 'ex74ftube.h'
3891: debug = .false.
3893: prat = p2/p1
3894: grat1 = (gamma - 1.0d+0) / (gamma + 1.0d+0)
3895: grat2 = (2.0d+0 * gamma) / (gamma + 1.0d+0)
3896: arat = a1/gamma
3898: xnum = sqrt ( grat2 / (prat + grat1) )
3900: v2 = arat * (prat - 1.0d+0) * xnum
3902: if (debug) then
3903: write(*,*)
3904: write(*,*) 'v2 = ',v2
3905: endif
3907: return
3908: end
3909: subroutine val3
3911: implicit none
3913: double precision prat, rpow, epow, p3t
3916: logical debug
3918: include 'ex74ftube.h'
3920: debug = .false.
3923: p3 = p2
3925: prat = p3/p4
3927: rpow = 1.0d+0 / gamma
3929: r3 = r4 * ( prat ** rpow )
3931: epow = (gamma - 1.0d+0) / gamma
3933: e3 = e4 * ( (p3/p4) ** epow )
3935: p3t = (gamma - 1.0d+0) * r3 * e3
3937: a3 = sqrt(gamma*p3/r3)
3939: if (debug) then
3940: write(*,*)
3941: write(*,*) 'a3 = ',a3
3942: write(*,*) 'r3 = ',r3
3943: write(*,*) 'e3 = ',e3
3944: write(*,*) 'p3 = ',p3
3945: write(*,*) 'p3t = ',p3t,' error = ',p3-p3t
3946: write(*,*)
3947: endif
3949: return
3950: end
3951: subroutine wval
3953: implicit none
3955: double precision prat, grat, xnum
3958: logical debug
3960: include 'ex74ftube.h'
3962: debug = .false.
3964: prat = p2/p1
3965: grat = (gamma + 1.0d+0) / (2.0d+0 * gamma)
3967: xnum = ( grat * (prat - 1.0d+0) ) + 1.0d+0
3968:
3969: sspd = a1 * sqrt(xnum)
3970:
3973: if (debug) then
3974: write(*,*)
3975: write(*,*) 'sspd = ',sspd
3976: endif
3978: return
3979: end