Actual source code: ex33f.F
1: >
2: program radhyd
3: ! "$Id: ex4f.F,v 1.39 1999/03/10 19:29:25 Vince Mousseau $";
4: !
5: ! Description: This example solves a nonlinear system on 1 processor with SNES.
6: ! We solve the Euler equations in one dimension.
7: ! The command line options include:
8: ! -dt <dt>, where <dt> indicates time step
9: ! -mx <xg>, where <xg> = number of grid points in the x-direction
10: ! -nstep <nstep>, where <nstep> = number of time steps
11: ! -debug <ndb>, where <ndb> = 0) no debug 1) debug
12: ! -pcnew <npc>, where <npc> = 0) no preconditioning 1) rad preconditioning
13: ! -probnum <probnum>, where <probnum> = 1) cyclic Riesner 2) dam break
14: ! -ihod <ihod>, where <ihod> = 1) upwind 2) quick 3) godunov
15: ! -ientro <ientro>, where <ientro> = 0) basic 1) entropy fix 2) hlle
16: ! -theta <theta>, where <theta> = 0-1 0-explicit 1-implicit
17: ! -hnaught <hnaught>, where <hnaught> = height of left side
18: ! -hlo <hlo>, where <hlo> = hieght of right side
19: ! -ngraph <ngraph>, where <ngraph> = number of time steps between graphics
20: ! -damfac <damfac>, where <damfac> = fractional downward change in hight
21: ! -dampit <ndamp>, where <ndamp> = 1 turn damping on
22: ! -gorder <gorder>, where <gorder> = spatial oerder of godunov
23: !
24: !
25: !
26: --------------------------------------------------------------------------
27: !
28: ! Shock tube example
29: !
30: ! In this example the application context is a Fortran integer array:
31: ! ctx(1) = shell preconditioner pressure matrix contex
32: ! ctx(2) = semi implicit pressure matrix
33: ! ctx(4) = xold - old time values need for time advancement
34: ! ctx(5) = mx - number of control volumes
35: ! ctx(6) = N - total number of unknowns
36: !
37: !
38: --------------------------------------------------------------------------
40: implicit none
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Include files
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! The following include statements are generally used in SNES Fortran
47: ! programs:
48: ! petsc.h - base PETSc routines
49: ! vec.h - vectors
50: ! mat.h - matrices
51: ! ksp.h - Krylov subspace methods
52: ! pc.h - preconditioners
53: ! snes.h - SNES interface
54: ! In addition, we need the following for use of PETSc drawing routines
55: ! draw.h - drawing routines
56: ! Other include statements may be needed if using additional PETSc
57: ! routines in a Fortran program, e.g.,
58: ! viewer.h - viewers
59: ! is.h - index sets
60: !
61: #include include/finclude/petsc.h
62: #include include/finclude/petscvec.h
63: #include include/finclude/petscis.h
64: #include include/finclude/petscdraw.h
65: #include include/finclude/petscmat.h
66: #include include/finclude/petscksp.h
67: #include include/finclude/petscpc.h
68: #include include/finclude/petscsnes.h
69: #include include/finclude/petscviewer.h
71: #include "comd.h"
72: #include "tube.h"
74: !
75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: ! Variable declarations
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: !
79: ! Variables:
80: ! snes - nonlinear solver
81: ! x,r - solution, residual vectors
82: ! J - Jacobian matrix
83: ! its - iterations for convergence
84: ! matrix_free - flag - 1 indicates matrix-free version
85: ! dt - time step size
86: ! draw - drawing context
87: !
88: PetscFortranAddr ctx(6)
89: integer nx,ny
90: SNES snes
91: KSP ksp
92: PC pc
93: Vec x,r
94: PetscViewer view0,view1,view2,
95: 1 view3, view4
96: Mat Psemi
97: integer matrix_free, flg, N, ierr, ngraph
98: integer nstep, ndt, size, rank, i
99: integer its, lits, totits, totlits
100: integer ndb, npc, ndamp, nwilson, ndtcon
101: double precision plotim
102: ! logical pcnew
104: double precision krtol,katol,kdtol
105: double precision natol,nrtol,nstol
106: integer kmit,nmit,nmf
109: ! Note: Any user-defined Fortran routines (such as FormJacobian)
110: ! MUST be declared as external.
112: external FormFunction, FormInitialGuess,FormDt,
113: & PCRadSetUp, PCRadApply, FormGraph,FormDampit
114: double precision eos
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: ! Initialize program
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: open (unit=87,file='Dt.out',status='unknown')
122: c
123: c start PETSc
124: c
125: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
126: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
127: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
129: if (size .ne. 1) then
130: if (rank .eq. 0) then
131: write(6,*) 'This is a uniprocessor example only!'
132: endif
133: SETERRQ(1,' ',ierr)
134: endif
136: ! Initialize problem parameters
138: debug = .false.
139: dampit = .false.
140: wilson = .true.
141: dtcon = .true.
142: pcnew = .true.
143: dtmax = 1.0d+2
144: dtmin = 1.00d-12
145: dt = 1.0d-2
146: mx = 100
147: nstep = 50
148: matrix_free = 1
149: probnum = 1
150: gorder = 1
152: tfinal = 1.0d+0
153: tplot = 0.2d+0
154: dtgrow = 1.05d+0
155: tcscal = 0.5d+0
156: hcscal = 0.5d+0
158: ihod = 3
159: ientro = 1
160: theta = 1.00d+0
161: pi = 3.14159d+0
163: zero = 0.0
164: ngraph = 10
166: ndb = 0
167: npc = 1
169: damfac = 0.9d+0
171: gamma = 1.25d+0
172: csubv = 1.0d+0 / (gamma - 1.0d+0)
174: v1 = 0.0d+0
175: v4 = 0.0d+0
177: e1 = 1.0d+0
178: e4 = 1.0d+0
180: r1 = 1.0d+0
181: r4 = 2.0d+0
183: ru1 = r1 * v1
184: ru4 = r4 * v4
186: et1 = r1 * ( (0.5d+0 * v1 * v1) + e1 )
187: et4 = r4 * ( (0.5d+0 * v4 * v4) + e4 )
189: p1 = eos(r1,ru1,et1)
190: p4 = eos(r4,ru4,et4)
192: a1 = sqrt(gamma*p1/r1)
193: a4 = sqrt(gamma*p4/r4)
195: erg0 = 1.0d+2
196: kappa0 = 1.0d+0
197: kappaa = -2.0d+0
198: kappab = 13.0d+0 / 2.0d+0
199: CVAM kappab = 2.5d+0
201: c
202: c load the command line options
203: c
204: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
205: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
206: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-nstep',nstep,flg,
207: & ierr)
208: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-debug',ndb,flg,
209: & ierr)
210: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-pcnew',npc,flg,
211: & ierr)
212: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ihod',ihod,flg,
213: & ierr)
214: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ientro',ientro,flg,
215: & ierr)
216: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',
217: & theta,flg,ierr)
218: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-ngraph',ngraph,flg,
219: & ierr)
220: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
221: & '-damfac',damfac,flg,ierr)
222: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dampit',ndamp,flg,
223: & ierr)
224: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
225: & '-wilson',nwilson,flg,ierr)
226: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-gorder',gorder,flg,
227: & ierr)
228: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
229: & '-probnum',probnum,flg,ierr)
230: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
231: & '-kappa0',kappa0,flg,ierr)
232: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
233: & '-erg0',erg0,flg,ierr)
234: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-dtcon',ndtcon,flg,
235: & ierr)
236: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
237: & '-tfinal',tfinal,flg,ierr)
238: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
239: & '-tplot',tplot,flg,ierr)
240: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
241: & '-dtgrow',dtgrow,flg,ierr)
242: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
243: & '-tcscal',tcscal,flg,ierr)
244: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
245: & '-hcscal',hcscal,flg,ierr)
246: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,
247: & '-dtmax',dtmax,flg,ierr)
249: if (ndamp .eq. 1) then
250: dampit = .true.
251: endif
253: if (nwilson .eq. 0) then
254: wilson = .false.
255: endif
257: if (ndb .eq. 1) then
258: debug = .true.
259: endif
261: if (npc .eq. 0) then
262: pcnew = .false.
263: endif
265: if (ndtcon .eq. 0) then
266: dtcon = .false.
267: endif
269: CVAM if (dt .ge. dtmax .or. dt .le. dtmin) then
270: CVAM if (rank .eq. 0) write(6,*) 'DT is out of range'
271: CVAM SETERRA(1,0,' ')
272: CVAM endif
274: N = mx*neq
276: ctx(5) = mx
277: ctx(6) = N
279: if (debug) then
280: write(*,*) 'mx = ',mx
281: endif
285: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: ! Create nonlinear solver context
287: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
291: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
292: ! Create vector data structures; set function evaluation routine
293: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295: CALL MatCreate(PETSC_COMM_WORLD, ctx(2), ierr)
296: CALL MatSetSizes(ctx(2),PETSC_DECIDE,PETSC_DECIDE,mx,mx,ierr)
297: CALL MatSetFromOptions(ctx(2),ierr)
298: if (debug) then
299: call MatGetSize(ctx(2),nx,ny,ierr)
300: write(*,*) 'number of rows = ',nx,' number of col = ',ny
301: endif
302: c
303: c full size vectors
304: c
305: CALL VecCreate(PETSC_COMM_WORLD,x,ierr)
306: CALL VecSetSizes(x,PETSC_DECIDE,N,ierr)
307: CALL VecSetType(x,VECMPI,ierr)
308: call VecSetFromOptions(x,ierr)
309: call VecDuplicate(x,r,ierr)
310: call VecDuplicate(x,ctx(4),ierr)
311: c
312: c set grid
313: c
314: dx = 2.0d+0/dfloat(mx)
315: xl0 = -1.0d+0 -(0.5d+0 * dx)
317: if (debug) then
318: write(*,*) 'dx = ',dx
319: endif
320:
322: ! Set function evaluation routine and vector. Whenever the nonlinear
323: ! solver needs to evaluate the nonlinear function, it will call this
324: ! routine.
325: ! - Note that the final routine argument is the user-defined
326: ! context that provides application-specific data for the
327: ! function evaluation routine.
329: call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
331: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332: ! Customize nonlinear solver; set runtime options
333: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type
336: <type>)
338: call SNESSetFromOptions(snes,ierr)
339: c
340: c set the line search function to damp the newton update.
341: c
342: ! if (dampit) then
343: ! call SNESSetLineSearch(snes,FormDampit,ctx,ierr)
344: ! endif
345: c
346: c get the linear solver info from the nonlinear solver
347: c
349: call SNESGetKSP(snes,ksp,ierr)
350: call KSPGetPC(ksp,pc,ierr)
351: ! call KSPGetKSP(ksp,ksp1,ierr)
352: CVAM call KSPSetType(ksp,KSPPREONLY,ierr)
353: call KSPSetType(ksp,KSPGMRES,ierr)
355: call KSPGetTolerances(ksp,krtol,katol,kdtol,kmit,ierr)
356: call SNESGetTolerances(snes,natol,nrtol,nstol,nmit,nmf,ierr)
358: write(*,*)
359: write(*,*)
360: write(*,*) 'Linear solver'
361: write(*,*)
362: write(*,*) 'rtol = ',krtol
363: write(*,*) 'atol = ',katol
364: write(*,*) 'dtol = ',kdtol
365: write(*,*) 'maxits = ',kmit
366: write(*,*)
367: write(*,*)
368: write(*,*) 'Nonlinear solver'
369: write(*,*)
370: write(*,*) 'rtol = ',nrtol
371: write(*,*) 'atol = ',natol
372: write(*,*) 'stol = ',nstol
373: write(*,*) 'maxits = ',nmit
374: write(*,*) 'max func = ',nmf
375: write(*,*)
376: write(*,*)
378: c
379: c Build shell based preconditioner if flag set
380: c
381: if (pcnew) then
382: call PCSetType(pc,PCSHELL,ierr)
383: call PCShellSetContext(pc,ctx,ierr)
384: call PCShellSetSetUp(pc,PCRadSetUp,ierr)
385: call PCShellSetApply(pc,PCRadApply,ierr)
386: endif
388: call PCCreate(PETSC_COMM_WORLD,ctx(1),ierr)
390: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391: ! Evaluate initial guess; then solve nonlinear system.
392: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
393: c
394: c initial counters
395: c
396: time = 0.0d+0
397: plotim = 0.0d+0
398: totits = 0
399: totlits = 0
401: ! Note: The user should initialize the vector, x, with the initial guess
402: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
403: ! to employ an initial guess of zero, the user should explicitly set
404: ! this vector to zero by calling VecSet().
406: call FormInitialGuess(x,ierr)
407: c
408: c open a window to plot results
409: c
410: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
411: & 'density',0,0,300,300,view0,ierr)
412: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
413: & 'velocity',320,0,300,300,view1,ierr)
414: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
415: & 'total energy',640,0,300,300,view2,ierr)
416: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
417: & 'temperature',0,360,300,300,view3,ierr)
418: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,
419: & 'pressure',320,360,300,300,view4,ierr)
420: c
421: c graph initial conditions
422: c
423: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
424: c
425: c copy x into xold
426: c
427: call VecCopy(x,ctx(4),ierr)
428: call FormDt(snes,x,ctx,ierr)
429: c
430: c################################
431: c
432: c TIME STEP LOOP BEGIN
433: c
434: c################################
435: c
436: ndt = 0
438: 10 if ( (ndt .le. nstep) .and. ((time + 1.0d-10) .lt. tfinal) ) then
440: if (debug) then
441: write(*,*)
442: write(*,*) 'start of time loop'
443: write(*,*)
444: write(*,*) 'ndt = ',ndt
445: write(*,*) 'nstep = ',nstep
446: write(*,*) 'time = ',time
447: write(*,*) 'tfinal = ',tfinal
448: write(*,*)
449: endif
451: ndt = ndt + 1
452: c
453: c increment time
454: c
455: time = time + dt
456: plotim = plotim + dt
457: c
458: c call the nonlinear solver
459: c
460: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
461: CALL SNESGetIterationNumber(snes,its,ierr)
462: c
463: c get the number of linear iterations used by the nonlinear solver
464: c
465: call SNESGetNumberLinearIterations(snes,lits,ierr)
467: if (debug) then
468: write(*,*) 'in radhyd ',ndt,'x'
469: call VecView(x,PETSC_VIEWER_STDOUT_SELF,ierr)
470: endif
471: c
472: c increment the counters
473: c
474: totits = totits + its
475: totlits = totlits + lits
476: c
477: c Compute new time step
478: c
479: call FormDt(snes,x,ctx,ierr)
480: c
481: c Draw contour plot of solution
482: c
483: if ( (mod(ndt,ngraph) .eq. 0) .or. (plotim .gt. tplot) )then
484:
485: plotim = plotim - tplot
488: if (rank .eq. 0) then
489: write(6,100) totits,totlits,ndt,dt,time
490: endif
491: 100 format('Newt = ',i7,' lin =',i7,' step =',i7,
492: & ' dt = ',e8.3,' time = ',e10.4)
493: c
494: c graph state conditions
495: c
496: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
498: endif
499: c
500: c copy x into xold
501: c
502: call VecCopy(x,ctx(4),ierr)
505: goto 10
507: endif
508: c
509: c################################
510: c
511: c TIME STEP LOOP END
512: c
513: c################################
514: c
516: c
517: c graph final conditions
518: c
519: call FormGraph(x,view0,view1,view2,view3,view4,ierr)
522: write(*,*)
523: write(*,*)
524: write(*,*) 'total Newton iterations = ',totits
525: write(*,*) 'total linear iterations = ',totlits
526: write(*,*) 'Average Newton per time step = ',
527: & dble(totits)/dble(ndt)
528: write(*,*) 'Average Krylov per Newton = ',
529: & dble(totlits)/dble(totits)
530: write(*,*)
531: write(*,*)
533: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
534: ! Free work space. All PETSc objects should be destroyed when they
535: ! are no longer needed.
536: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
539: call MatDestroy(ctx(2),ierr)
540: call VecDestroy(x,ierr)
541: call VecDestroy(ctx(4),ierr)
542: call VecDestroy(r,ierr)
543: call SNESDestroy(snes,ierr)
544: call PetscViewerDestroy(view0,ierr)
545: call PetscViewerDestroy(view1,ierr)
546: call PetscViewerDestroy(view2,ierr)
547: call PetscViewerDestroy(view3,ierr)
548: call PetscViewerDestroy(view4,ierr)
550: call PCDestroy(ctx(1),ierr)
552: call PetscFinalize(ierr)
554: close(87)
556: stop
557: end
558: subroutine ApplicationDampit(x,deltx,w,ierr)
559: ! ---------------------------------------------------------------------
560: !
561: ! ApplicationDampit - Damps the newton update, called by
562: ! the higher level routine FormDampit().
563: !
564: ! Input Parameter:
565: ! x - current iterate
566: ! deltx - update
567: !
568: ! Output Parameters:
569: ! x - new iterate
570: ! ierr - error code
571: !
572: ! Notes:
573: ! This routine only damps the density. May want to add energy
574: ! in the future
575: !
577: implicit none
579: ! Common blocks:
580: #include "comd.h"
582: ! Input/output variables:
583: PetscScalar x(mx*neq),deltx(mx*neq),
584: 1 w(mx*neq)
585: integer ierr
587: ! Local variables:
588: double precision facmin, fac, newx, xmin, se, dse
589: double precision u,en,rn,run
590: integer i, jr, jru, je
592: 0
594: if (debug) then
595: write(*,*) 'begin damping'
596: do i = 1,mx*neq
597: write(*,*)'x(',i,') = ',x(i)
598: enddo
599: write(*,*)
600: do i = 1,mx*neq
601: write(*,*)'deltx(',i,') = ',deltx(i)
602: enddo
603: endif
605: facmin = 1.0d+0
606: c
607: c set the scale factor
608: c
609: do i=1,mx
610: c
611: c set pointers
612: c
613: jr = (neq*i) - 2
614: jru = (neq*i) - 1
615: je = (neq*i)
616: c
617: c make sure dencity stayes positive
618: c
619: newx = x(jr) - deltx(jr)
620: xmin = damfac * x(jr)
622: if (newx .lt. xmin) then
623: fac = (1.0d+0 - damfac)*x(jr)/deltx(jr)
624: if (fac .lt. facmin) then
625: if (debug) then
626: write(*,*) 'density', i, damfac,facmin,fac,x(jr),deltx(jr)
627: endif
628: facmin = fac
629: endif
630: endif
631: c
632: c make sure Total energy stayes positive
633: c
634: newx = x(je) - deltx(je)
635: xmin = damfac * x(je)
637: if (newx .lt. xmin) then
638: fac = (1.0d+0 - damfac)*x(je)/deltx(je)
639: if (fac .lt. facmin) then
640: if (debug) then
641: write(*,*) 'energy T',i, damfac,facmin,fac,x(je),deltx(je)
642: endif
643: facmin = fac
644: endif
645: endif
646: c
647: c make sure specific internal energy stayes positive
648: c
649:
650: u = x(jru)/x(jr)
651: se = (x(je)/x(jr)) - (0.5d+0 * u * u)
653: en = x(je) - deltx(je)
654: rn = x(jr) - deltx(jr)
655: run = x(jru) - deltx(jru)
657: dse = se - ( (en/rn) - (0.5d+0 * (run/rn) * (run/rn)) )
660: newx = se - dse
661: xmin = damfac * se
663: if (newx .lt. xmin) then
664: fac = (1.0d+0 - damfac) * se / dse
665: if (fac .lt. facmin) then
666: if (debug) then
667: write(*,*) 'se',i, damfac,facmin,fac,se,dse
668: endif
669: facmin = fac
670: endif
671: endif
673: enddo
674: c
675: c write out warning
676: c
677: if (facmin .lt. 1.0d+0) then
678: write(*,*) 'facmin = ',facmin, damfac,time
679: c
680: c scale the vector
681: c
682: do i=1,neq*mx
683: w(i) = x(i) - (facmin * deltx(i))
684: enddo
685: else
686: do i=1,neq*mx
687: w(i) = x(i) - deltx(i)
688: enddo
689: endif
691: if (debug) then
692: write(*,*) 'end damping'
693: do i = 1,mx*neq
694: write(*,*) 'w(',i,') = ',w(i)
695: enddo
696: endif
698: return
699: end
700: subroutine ApplicationDt(x,xold,ierr)
701: ! ---------------------------------------------------------------------
702: !
703: ! ApplicationDt - compute CFL numbers. Called by
704: ! the higher level routine FormDt().
705: !
706: ! Input Parameter:
707: ! x - local vector data
708: !
709: ! Output Parameters:
710: ! ierr - error code
711: !
712: ! Notes:
713: ! This routine uses standard Fortran-style computations over a 2-dim
714: array.
715: !
717: implicit none
719: ! Common blocks:
720: #include "comd.h"
721: #include "tube.h"
723: ! Input/output variables:
724: PetscScalar x(mx*neq), xold(mx*neq)
725: integer ierr
727: ! Local variables:
728: integer i, jr, jru, je
729: c
730: c new
731: c
732: double precision rhoee, rhoe, rhop, rhow, rhoww,
733: & rhouee, rhoue, rhoup, rhouw, rhouww,
734: & ergee, erge, ergp, ergw, ergww,
735: & vele, velp, velw
736: double precision pressp,sndp, vrad, vradn, vradd, tcfl, hcfl
737: double precision tcflg, hcflg, dtt, dth
738: double precision te, tp, tw
739: double precision ue, up, uw
740: double precision see, sep, sew
741: c
742: c old
743: c
744: double precision rhooee, rhooe, rhoop, rhoow, rhooww
745: double precision rhouoee, rhouoe, rhouop, rhouow, rhouoww
746: double precision ergoee, ergoe, ergop, ergow, ergoww
747: double precision veloe, velop, velow
748: double precision uop, seop, top
749: double precision dtold, dttype
750: c
751: c functions
752: c
753: double precision eos
755: dtold = dt
757: 0
759: if (debug) then
760: write(*,*) 'in start dt'
761: do i = 1,mx*neq
762: write(*,*)'x(',i,') = ',x(i)
763: enddo
764: write(*,*) 'tfinal = ',tfinal
765: write(*,*) 'time = ',time
766: write(*,*) 'dt = ',dt
767: write(*,*) 'dtmax = ',dtmax
768: endif
770: sndp = -1.0d+20
771: vradn = 0.0d+0
772: vradd = 0.0d+0
774: c
775: c################################
776: c
777: c loop over all cells begin
778: c
779: c################################
780: c
781: do i=1,mx
782: c
783: c set pointers
784: c
785: jr = (neq*i) - 2
786: jru = (neq*i) - 1
787: je = (neq*i)
788: c
789: c
790: c set scalars
791: c
792: call Setpbc(i,x,
793: & rhoee, rhoe, rhop, rhow, rhoww,
794: & rhouee, rhoue, rhoup, rhouw, rhouww,
795: & ergee, erge, ergp, ergw, ergww,
796: & vele, velp, velw)
797: c
798: c compute temperatures
799: c
800: uw = rhouw / rhow
801: up = rhoup / rhop
802: ue = rhoue / rhoe
804: see = (erge/rhoe) - (0.5d+0 * ue * ue)
805: sep = (ergp/rhop) - (0.5d+0 * up * up)
806: sew = (ergw/rhow) - (0.5d+0 * uw * uw)
808: te = see / csubv
809: tp = sep / csubv
810: tw = sew / csubv
811: c
812: c compute old temperature
813: c
815: call Setpbc(i,xold,
816: & rhooee, rhooe, rhoop, rhoow, rhooww,
817: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
818: & ergoee, ergoe, ergop, ergow, ergoww,
819: & veloe, velop, velow)
821: call Setpbcn(rhooee, rhooe, rhoop, rhoow, rhooww,
822: & rhouoee, rhouoe, rhouop, rhouow, rhouoww,
823: & ergoee, ergoe, ergop, ergow, ergoww,
824: & veloe, velop, velow, i)
826: uop = rhouop / rhoop
828: seop = (ergop/rhoop) - (0.5d+0 * uop * uop)
830: top = seop / csubv
831: c
832: c compute thermal cfl
833: c
834: vradn = vradn + (abs(tp - top)/dt)
835: vradd = vradd + (abs(te - tw) / (2.0d+0 * dx) )
836: c
837: c compute hydro cfl
838: c
840: pressp = eos(rhop, rhoup, ergp)
841: sndp = max(sndp,sqrt( (gamma*pressp) / rhop ))
843: enddo
844: c
845: c################################
846: c
847: c loop over all cells end
848: c
849: c################################
850: c
852: vrad = vradn / vradd
854: tcfl = (vrad * dt) / dx
855: hcfl = (sndp * dt) / dx
857: dtt = max(dx/vrad,1.0d-7)
858: dtt = tcscal * dtt
860: dth = hcscal * dx / sndp
862: if (.not. dtcon) then
863: dt = min (dth,dtt,dt*dtgrow)
864: if (dt .lt. dtmin) then
865: dt = dtmin
866: endif
867: if (dt .gt. dtmax) then
868: dt = dtmax
869: endif
870: if ( (time + dt) .gt. tfinal) then
871: dt = tfinal - time
872: endif
874: if (dt .eq. dth) then
875: dttype = 1.0d+0
876: elseif (dt .eq. dtt) then
877: dttype = 2.0d+0
878: elseif (dt .eq. dtold*dtgrow) then
879: dttype = 3.0d+0
880: elseif (dt .eq. dtmax) then
881: dttype = 4.0d+0
882: elseif (dt .eq. dtmin) then
883: dttype = 5.0d+0
884: elseif (dt .eq. tfinal - time) then
885: dttype = 6.0
886: else
887: dttype = -1.0d+0
888: endif
890: endif
891:
892:
893: write(87,1000) time,dt,dth/hcscal,dtt/tcscal
894: write(88,1000) time,dttype
896: 1000 format(4(2x,e18.9))
898: if (debug) then
899: write(*,*) 'thermal cfl = ',tcfl,'hydro cfl = ',hcfl
900: write(*,*) 'dtt = ',dtt,' dth = ',dth
901: write(*,*) 'tfinal = ',tfinal
902: write(*,*) 'time = ',time
903: write(*,*) 'dt = ',dt
904: write(*,*) 'dtmax = ',dtmax
905: write(*,*)
906: endif
908: return
909: end
910: subroutine ApplicationExact(x,ierr)
911: ! ---------------------------------------------------------------------
912: !
913: ! ApplicationExact - Computes exact solution, called by
914: ! the higher level routine FormExact().
915: !
916: ! Input Parameter:
917: ! x - local vector data
918: !
919: ! Output Parameters:
920: ! x - initial conditions
921: ! ierr - error code
922: !
923: ! Notes:
924: ! This routine uses standard Fortran-style computations over a 1-dim
925: array.
926: !
928: implicit none
930: ! Common blocks:
932: #include "comd.h"
934: ! Input/output variables:
935: PetscScalar x(mx)
936: integer ierr
938: ! Local variables:
939: integer i
940: double precision xloc
941: PetscScalar rexact
944: ! Set parameters
946: 0
948: do i = 1,mx
950: xloc = xl0 + (dble(i) * dx)
951: x(i) = rexact(xloc,time)
953: enddo
955: return
956: end
957: subroutine ApplicationFunction(x,f,xold,ierr)
958: ! ---------------------------------------------------------------------
959: !
960: ! ApplicationFunction - Computes nonlinear function, called by
961: ! the higher level routine FormFunction().
962: !
963: ! Input Parameter:
964: ! x - local vector data
965: !
966: ! Output Parameters:
967: ! f - local vector data, f(x)
968: ! ierr - error code
969: !
970: ! Notes:
971: ! This routine uses standard Fortran-style computations over a 2-dim
972: array.
973: !
975: implicit none
977: ! Common blocks:
978: #include "comd.h"
980: ! Input/output variables:
981: PetscScalar x(mx*neq),f(mx*neq),
982: 1 xold(mx*neq)
983: integer ierr
985: ! Local variables:
986: integer i, jr, jru, je
987: double precision rhoee, rhoe, rhop, rhow, rhoww,
988: & rhouee, rhoue, rhoup, rhouw, rhouww,
989: & ergee, erge, ergp, ergw, ergww,
990: & vele, velp, velw
992: double precision cont, energy, mom
994: 0
996: if (debug) then
997: write(*,*) 'in function'
998: do i = 1,mx*neq
999: write(*,*)'x(',i,') = ',x(i)
1000: enddo
1001: endif
1002: c
1003: c################################
1004: c
1005: c loop over all cells begin
1006: c
1007: c################################
1008: c
1009: do i=1,mx
1010: c
1011: c set pointers
1012: c
1013: jr = (neq*i) - 2
1014: jru = (neq*i) - 1
1015: je = (neq*i)
1016: c
1017: c
1018: c set scalars
1019: c
1020: call Setpbc(i,x,
1021: & rhoee, rhoe, rhop, rhow, rhoww,
1022: & rhouee, rhoue, rhoup, rhouw, rhouww,
1023: & ergee, erge, ergp, ergw, ergww,
1024: & vele, velp, velw)
1025: c
1026: c compute functions
1027: c
1029: f(jr) = cont(rhoee, rhoe, rhop, rhow, rhoww,
1030: & rhouee, rhoue, rhoup, rhouw, rhouww,
1031: & ergee, erge, ergp, ergw, ergww,
1032: & i,xold)
1035: f(jru) = mom(rhoee, rhoe, rhop, rhow, rhoww,
1036: & rhouee, rhoue, rhoup, rhouw, rhouww,
1037: & ergee, erge, ergp, ergw, ergww,
1038: & i,xold)
1041: f(je) = energy(rhoee, rhoe, rhop, rhow, rhoww,
1042: & rhouee, rhoue, rhoup, rhouw, rhouww,
1043: & ergee, erge, ergp, ergw, ergww,
1044: & i,xold)
1046: if (debug) then
1047: write(*,*)
1048: write(*,*) i,jr,jru,je,'res,r,ru,e'
1049: write(*,*) f(jr),f(jru),f(je)
1050: write(*,*)
1051: endif
1053: enddo
1054: c
1055: c################################
1056: c
1057: c loop over all cells end
1058: c
1059: c################################
1060: c
1062: if (debug) then
1063: write(*,*) 'in function'
1064: do i = 1,mx*neq
1065: write(*,*) 'f(',i,') = ',f(i)
1066: enddo
1067: endif
1069: return
1070: end
1071: subroutine ApplicationInitialGuess(x,ierr)
1072: ! ---------------------------------------------------------------------
1073: !
1074: ! ApplicationInitialGuess - Computes initial approximation, called by
1075: ! the higher level routine FormInitialGuess().
1076: !
1077: ! Input Parameter:
1078: ! x - local vector data
1079: !
1080: ! Output Parameters:
1081: ! x - initial conditions
1082: ! ierr - error code
1083: !
1084: ! Notes:
1085: ! This routine uses standard Fortran-style computations over a 1-dim
1086: array.
1087: !
1089: implicit none
1091: ! Common blocks:
1093: #include "comd.h"
1094: #include "tube.h"
1096: ! Input/output variables:
1097: PetscScalar x(mx*neq)
1098: integer ierr
1100: ! Local variables:
1101: integer i, j, jr, jru, je
1102: double precision xloc, re, ee, ve
1103: double precision wid, efloor
1104: PetscScalar uexact, rexact, eexact
1107: CVAM efloor = max(1.0d-1,1.0d-3 * erg0)
1108: efloor = 1.0d-1
1109: CVAM wid = max(1.0d-2,dx)
1110: wid = 1.0d-2
1112: ! Set parameters
1114: 0
1116: do i = 1,mx
1118: jr = (neq*i) - 2
1119: jru = (neq*i) - 1
1120: je = (neq*i)
1122: xloc = xl0 + (dble(i) * dx)
1124: if (probnum .eq. 1) then
1125: re = rexact(xloc,zero)
1126: ve = uexact(xloc,zero)
1127: ee = eexact(xloc,zero)
1128: else
1129: re = 1.0d+0
1130: ve = 0.0d+0
1131: ee = efloor + (erg0 * exp(-(xloc*xloc)/(wid*wid)))
1132: endif
1134: x(jr) = re
1135: x(jru) = re * ve
1136: x(je) = re * ( (0.5d+0 * ve * ve) + ee )
1138: if (debug) then
1139: write(*,100) i,jr,jru,je,xloc,x(jr),x(jru),x(je)
1140: 100 format(i3,2x,i3,2x,i3,2x,i3,4(2x,e12.5))
1141: endif
1143: enddo
1145: call exact0
1146: call eval2
1147: call rval2
1148: call wval
1149: call uval2
1150: v3 = v2
1151: call val3
1153: a1 = sqrt(gamma*p1/r1)
1154: a2 = sqrt(gamma*p2/r2)
1155: a3 = sqrt(gamma*p3/r3)
1156: a4 = sqrt(gamma*p4/r4)
1158: write(*,1000) r1,r2,r3,r4
1159: write(*,2000) p1,p2,p3,p4
1160: write(*,3000) e1,e2,e3,e4
1161: write(*,4000) a1,a2,a3,a4
1162: write(*,*)
1164: 1000 format ('rhos ',4(f12.6))
1165: 2000 format ('pressures ',4(f12.6))
1166: 3000 format ('energies ',4(f12.6))
1167: 4000 format ('sound ',4(f12.6))
1170: return
1171: end
1172: subroutine ApplicationXmgr(x,ivar,ierr)
1173: ! ---------------------------------------------------------------------
1174: !
1175: ! ApplicationXmgr - Sets the Xmgr output called from
1176: ! the higher level routine FormXmgr().
1177: !
1178: ! Input Parameter:
1179: ! x - local vector data
1180: !
1181: ! Output Parameters:
1182: ! x - initial conditions
1183: ! ierr - error code
1184: !
1185: ! Notes:
1186: ! This routine uses standard Fortran-style computations over a 1-dim
1187: array.
1188: !
1190: implicit none
1192: ! Common blocks:
1194: #include "comd.h"
1196: ! Input/output variables:
1197: PetscScalar x(mx)
1198: integer ivar,ierr
1200: ! Local variables:
1201: integer i
1202: double precision xloc, sum
1203: PetscScalar rexact
1204: integer iplotnum(5)
1205: save iplotnum
1206: character*8 grfile
1209: data iplotnum / -1,-1,-1,-1,-1 /
1213: ! Set parameters
1215: iplotnum(ivar) = iplotnum(ivar) + 1
1216: 0
1218: if (ivar .eq. 1) then
1219: write(grfile,4000) iplotnum(ivar)
1220: 4000 format('Xmgrr',i3.3)
1221: elseif (ivar .eq. 2) then
1222: write(grfile,5000) iplotnum(ivar)
1223: 5000 format('Xmgru',i3.3)
1224: elseif (ivar .eq. 3) then
1225: write(grfile,6000) iplotnum(ivar)
1226: 6000 format('Xmgre',i3.3)
1227: elseif (ivar .eq. 4) then
1228: write(grfile,7000) iplotnum(ivar)
1229: 7000 format('Xmgrt',i3.3)
1230: else
1231: write(grfile,8000) iplotnum(ivar)
1232: 8000 format('Xmgrp',i3.3)
1233: endif
1235: open (unit=44,file=grfile,status='unknown')
1237: do i = 1,mx
1239: xloc = xl0 + (dble(i) * dx)
1240: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1241: write(44,1000) xloc, x(i), rexact(xloc,time)
1242: else
1243: write(44,1000) xloc, x(i)
1244: endif
1246: enddo
1248: 1000 format(3(e18.12,2x))
1249: close(44)
1251: if ( (ivar .eq. 1) .and. (probnum .eq. 1) ) then
1252: sum = 0.0d+0
1253: do i = 1,mx
1254: xloc = xl0 + (dble(i) * dx)
1255: sum = sum + (x(i) - rexact(xloc,time)) ** 2
1256: enddo
1257: sum = sqrt(sum)
1259: write(*,*)
1260: write(*,*) 'l2norm of the density error is',sum
1261: write(*,*)
1262: endif
1265: return
1266: end
1267: subroutine FormDampit(snes,ctx,x,f,g,y,w,
1268: & fnorm,ynorm,gnorm,flag,ierr)
1269: ! ---------------------------------------------------------------------
1270: !
1271: ! FormDampit - damps the Newton update
1272: !
1273: ! Input Parameters:
1274: ! snes - the SNES context
1275: ! x - current iterate
1276: ! f - residual evaluated at x
1277: ! y - search direction (containes new iterate on output)
1278: ! w - work vector
1279: ! fnorm - 2-norm of f
1280: !
1281: ! In this example the application context is a Fortran integer array:
1282: ! ctx(1) = shell preconditioner pressure matrix contex
1283: ! ctx(2) = semi implicit pressure matrix
1284: ! ctx(4) = xold - old time values need for time advancement
1285: ! ctx(5) = mx - number of control volumes
1286: ! ctx(6) = N - total number of unknowns
1287: !
1288: ! Output Parameter:
1289: ! g - residual evaluated at new iterate y
1290: ! y - new iterate (contains search direction on input
1291: ! gnorm - 2-norm of g
1292: ! ynorm - 2-norm of search length
1293: ! flag - set to 0 if the line search succeeds; -1 on failure
1294: !
1295: ! Notes:
1296: ! This routine serves as a wrapper for the lower-level routine
1297: ! "ApplicationDampit", where the actual computations are
1298: ! done using the standard Fortran style of treating the local
1299: ! vector data as a multidimensional array over the local mesh.
1300: ! This routine merely accesses the local vector data via
1301: ! VecGetArray() and VecRestoreArray().
1302: !
1303: implicit none
1305: #include include/finclude/petsc.h
1306: #include include/finclude/petscvec.h
1307: #include include/finclude/petscsnes.h
1309: ! Input/output variables:
1310: SNES snes
1311: Vec x, f, g, y, w
1312: PetscFortranAddr ctx(*)
1313: PetscScalar fnorm, ynorm, gnorm
1314: integer ierr, flag
1316: ! Common blocks:
1318: #include "comd.h"
1320: ! Local variables:
1322: ! Declarations for use with local arrays:
1323: PetscScalar lx_v(0:1),ly_v(0:1),
1324: 1 lw_v(0:1)
1325: PetscOffset lx_i,ly_i,lw_i
1327: c
1328: c set ynorm
1329: c
1330: call VecNorm(y,NORM_2,ynorm,ierr)
1331: c
1332: c copy x to w
1333: c
1334: call VecCopy(x,w,ierr)
1335: c
1336: c get pointers to x, y, w
1337: c
1339: call VecGetArray(x,lx_v,lx_i,ierr)
1340: call VecGetArray(y,ly_v,ly_i,ierr)
1341: call VecGetArray(w,lw_v,lw_i,ierr)
1342: c
1343: c Compute Damping
1344: c
1345: call ApplicationDampit(lx_v(lx_i),ly_v(ly_i),lw_v(lw_i),ierr)
1346: c
1347: c Restore vectors x, y, w
1348: c
1349: call VecRestoreArray(x,lx_v,lx_i,ierr)
1350: call VecRestoreArray(y,ly_v,ly_i,ierr)
1351: call VecRestoreArray(w,lw_v,lw_i,ierr)
1352: c
1353: c copy w to y
1354: c
1355: call VecCopy(w,y,ierr)
1356: c
1357: c compute new residual
1358: c
1359: call SNESComputeFunction(snes,y,g,ierr)
1360: call VecNorm(g,NORM_2,gnorm,ierr)
1361: flag = 0
1363: if (debug) then
1364: write(*,*) 'form damp ynorm = ',ynorm
1365: write(*,*) 'gnorm = ',gnorm
1366: endif
1368: return
1369: end
1370: subroutine FormDt(snes,x,ctx,ierr)
1371: ! ---------------------------------------------------------------------
1372: !
1373: ! FormDt - Compute CFL numbers
1374: !
1375: ! Input Parameters:
1376: ! snes - the SNES context
1377: ! x - input vector
1378: !
1379: ! In this example the application context is a Fortran integer array:
1380: ! ctx(1) = shell preconditioner pressure matrix contex
1381: ! ctx(2) = semi implicit pressure matrix
1382: ! ctx(4) = xold - old time values need for time advancement
1383: ! ctx(5) = mx - number of control volumes
1384: ! ctx(6) = N - total number of unknowns
1385: !
1386: !
1387: ! Notes:
1388: ! This routine serves as a wrapper for the lower-level routine
1389: ! "ApplicationDt", where the actual computations are
1390: ! done using the standard Fortran style of treating the local
1391: ! vector data as a multidimensional array over the local mesh.
1392: ! This routine merely accesses the local vector data via
1393: ! VecGetArray() and VecRestoreArray().
1394: !
1395: implicit none
1397: #include include/finclude/petsc.h
1398: #include include/finclude/petscvec.h
1399: #include include/finclude/petscsnes.h
1401: ! Input/output variables:
1402: SNES snes
1403: Vec x
1404: PetscFortranAddr ctx(*)
1405: integer ierr
1407: ! Common blocks:
1409: #include "comd.h"
1411: ! Local variables:
1413: ! Declarations for use with local arrays:
1414: PetscScalar lx_v(0:1)
1415: PetscOffset lx_i
1416: PetscScalar lxold_v(0:1)
1417: PetscOffset lxold_i
1419: c
1420: c get pointers to x, xold
1421: c
1423: call VecGetArray(x,lx_v,lx_i,ierr)
1424: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1425: c
1426: c Compute function
1427: c
1428: call ApplicationDt(lx_v(lx_i),lxold_v(lxold_i),ierr)
1429: c
1430: c Restore vectors x, xold
1431: c
1432: call VecRestoreArray(x,lx_v,lx_i,ierr)
1433: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1435: return
1436: end
1437: subroutine FormExact(x,ierr)
1438: ! ---------------------------------------------------------------------
1439: !
1440: ! FormExact - Forms exact solution
1441: !
1442: ! Input Parameter:
1443: ! x - vector
1444: !
1445: ! Output Parameters:
1446: ! x - vector
1447: ! ierr - error code
1448: !
1449: ! Notes:
1450: ! This routine serves as a wrapper for the lower-level routine
1451: ! "ApplicationExact", where the actual computations are
1452: ! done using the standard Fortran style of treating the local
1453: ! vector data as a multidimensional array over the local mesh.
1454: ! This routine merely accesses the local vector data via
1455: ! VecGetArray() and VecRestoreArray().
1456: !
1457: implicit none
1459: #include include/finclude/petsc.h
1460: #include include/finclude/petscvec.h
1461: #include include/finclude/petscmat.h
1462: #include include/finclude/petscsnes.h
1464: ! Input/output variables:
1465: Vec x
1466: integer ierr
1468: ! Declarations for use with local arrays:
1469: PetscScalar lx_v(0:1)
1470: PetscOffset lx_i
1472: 0
1474: c
1475: c get a pointer to x
1476: c
1477: call VecGetArray(x,lx_v,lx_i,ierr)
1478: c
1479: c Compute initial guess
1480: c
1481: call ApplicationExact(lx_v(lx_i),ierr)
1482: c
1483: c Restore vector x
1484: c
1485: call VecRestoreArray(x,lx_v,lx_i,ierr)
1487: return
1488: end
1489: subroutine FormFunction(snes,x,f,ctx,ierr)
1490: ! ---------------------------------------------------------------------
1491: !
1492: ! FormFunction - Evaluates nonlinear function, f(x).
1493: !
1494: ! Input Parameters:
1495: ! snes - the SNES context
1496: ! x - input vector
1497: !
1498: ! In this example the application context is a Fortran integer array:
1499: ! ctx(1) = shell preconditioner pressure matrix contex
1500: ! ctx(2) = semi implicit pressure matrix
1501: ! ctx(4) = xold - old time values need for time advancement
1502: ! ctx(5) = mx - number of control volumes
1503: ! ctx(6) = N - total number of unknowns
1504: !
1505: ! Output Parameter:
1506: ! f - vector with newly computed function
1507: !
1508: ! Notes:
1509: ! This routine serves as a wrapper for the lower-level routine
1510: ! "ApplicationFunction", where the actual computations are
1511: ! done using the standard Fortran style of treating the local
1512: ! vector data as a multidimensional array over the local mesh.
1513: ! This routine merely accesses the local vector data via
1514: ! VecGetArray() and VecRestoreArray().
1515: !
1516: implicit none
1518: #include include/finclude/petsc.h
1519: #include include/finclude/petscvec.h
1520: #include include/finclude/petscsnes.h
1522: ! Input/output variables:
1523: SNES snes
1524: Vec x, f
1525: PetscFortranAddr ctx(*)
1526: integer ierr
1528: ! Common blocks:
1530: #include "comd.h"
1532: ! Local variables:
1534: ! Declarations for use with local arrays:
1535: PetscScalar lx_v(0:1), lf_v(0:1)
1536: PetscOffset lx_i, lf_i
1537: PetscScalar lxold_v(0:1)
1538: PetscOffset lxold_i
1540: c
1541: c get pointers to x, f, xold
1542: c
1544: call VecGetArray(x,lx_v,lx_i,ierr)
1545: call VecGetArray(f,lf_v,lf_i,ierr)
1546: call VecGetArray(ctx(4),lxold_v,lxold_i,ierr)
1547: c
1548: c Compute function
1549: c
1550: call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),
1551: & lxold_v(lxold_i),ierr)
1552: c
1553: c Restore vectors x, f, xold
1554: c
1555: call VecRestoreArray(x,lx_v,lx_i,ierr)
1556: call VecRestoreArray(f,lf_v,lf_i,ierr)
1557: call VecRestoreArray(ctx(4),lxold_v,lxold_i,ierr)
1558: c
1559: c something to do with profiling
1560: c
1561: call PetscLogFlops(11*mx,ierr)
1563: return
1564: end
1565: subroutine FormGraph(x,view0,view1,view2,view3,view4,ierr)
1566: ! ---------------------------------------------------------------------
1567: !
1568: ! FormGraph - Forms Graphics output
1569: !
1570: ! Input Parameter:
1571: ! x - vector
1572: ! time - scalar
1573: !
1574: ! Output Parameters:
1575: ! ierr - error code
1576: !
1577: ! Notes:
1578: ! This routine serves as a wrapper for the lower-level routine
1579: ! "ApplicationXmgr", where the actual computations are
1580: ! done using the standard Fortran style of treating the local
1581: ! vector data as a multidimensional array over the local mesh.
1582: ! This routine merely accesses the local vector data via
1583: ! VecGetArray() and VecRestoreArray().
1584: !
1585: implicit none
1587: #include include/finclude/petsc.h
1588: #include include/finclude/petscvec.h
1589: #include include/finclude/petscmat.h
1590: #include include/finclude/petscsnes.h
1592: #include "comd.h"
1593: #include "tube.h"
1595: ! Input/output variables:
1596: Vec x
1597: integer ierr
1599: ! Declarations for use with local arrays:
1600: IS rfrom,rto,
1601: 1 rufrom, ruto, efrom, eto
1602: Vec rval
1603: Vec uval
1604: Vec ruval
1605: Vec eval
1606: Vec seval
1607: Vec pval
1608: Vec kval
1609: Vec tval
1610: Vec steval
1611: VecScatter scatter
1612: PetscViewer view0,view1,
1613: 1 view2, view3, view4
1614: double precision minus1, l2err, gm1, csubvi
1617: csubvi = 1.0d+0 / csubv
1618: gm1 = gamma - 1.0d+0
1619: 0
1620: minus1 = -1.0d+0
1621: c
1622: c graphics vectors
1623: c
1624: CALL VecCreate(PETSC_COMM_WORLD,rval,ierr)
1625: CALL VecSetSizes(rval,PETSC_DECIDE,mx,ierr)
1626: CALL VecSetType(rval,VECMPI,ierr)
1627: call VecSetFromOptions(rval,ierr)
1628: call VecDuplicate(rval,uval,ierr)
1629: call VecDuplicate(rval,ruval,ierr)
1630: call VecDuplicate(rval,eval,ierr)
1631: call VecDuplicate(rval,seval,ierr)
1632: call VecDuplicate(rval,pval,ierr)
1633: call VecDuplicate(rval,kval,ierr)
1634: call VecDuplicate(rval,tval,ierr)
1635: call VecDuplicate(rval,steval,ierr)
1636: c
1637: c create index sets for vector scatters
1638: c
1639: call ISCreateStride(PETSC_COMM_WORLD,mx,0,neq,rfrom, ierr)
1640: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, rto, ierr)
1641: call ISCreateStride(PETSC_COMM_WORLD,mx,1,neq,rufrom,ierr)
1642: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, ruto, ierr)
1643: call ISCreateStride(PETSC_COMM_WORLD,mx,2,neq,efrom, ierr)
1644: call ISCreateStride(PETSC_COMM_WORLD,mx,0,1, eto, ierr)
1646: c
1647: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1648: c
1649: c
1650: c load rval from x
1651: c
1652: call VecScatterCreate(x,rfrom,rval,rto,scatter,ierr)
1653: call VecScatterBegin(scatter,x,rval,INSERT_VALUES,
1654: & SCATTER_FORWARD,ierr)
1655: call VecScatterEnd(scatter,x,rval,INSERT_VALUES,
1656: & SCATTER_FORWARD,ierr)
1657: call VecScatterDestroy(scatter,ierr)
1658: c
1659: c plot rval vector
1660: c
1661: call VecView(rval,view0,ierr)
1662: c
1663: c make xmgr plot of rval
1664: c
1665: call FormXmgr(rval,1,ierr)
1666: c
1667: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1668: c
1669: c
1670: c load eval from x
1671: c
1672: call VecScatterCreate(x,efrom,eval,eto,scatter,ierr)
1673: call VecScatterBegin(scatter,x,eval,INSERT_VALUES,
1674: & SCATTER_FORWARD,ierr)
1675: call VecScatterEnd(scatter,x,eval,INSERT_VALUES,
1676: & SCATTER_FORWARD,ierr)
1677: call VecScatterDestroy(scatter,ierr)
1678: c
1679: c plot eval vector
1680: c
1681: call VecView(eval,view2,ierr)
1682: c
1683: c make xmgr plot of eval
1684: c
1685: call FormXmgr(eval,3,ierr)
1686: c
1687: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1688: c
1690: c
1691: c load ruval from x
1692: c
1693: call VecScatterCreate(x,rufrom,ruval,ruto,scatter,ierr)
1694: call VecScatterBegin(scatter,x,ruval,INSERT_VALUES,
1695: & SCATTER_FORWARD,ierr)
1696: call VecScatterEnd(scatter,x,ruval,INSERT_VALUES,
1697: & SCATTER_FORWARD,ierr)
1698: call VecScatterDestroy(scatter,ierr)
1699: c
1700: c create u = ru / r
1701: c
1702: call VecPointwiseDivide(ruval,rval,uval,ierr)
1703: c
1704: c plot uval vector
1705: c
1706: call VecView(uval,view1,ierr)
1707: c
1708: c make xmgr plot of uval
1709: c
1710: call FormXmgr(uval,2,ierr)
1712: c
1713: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1714: c
1716: call VecPointwiseMult(kval,uval,uval,ierr)
1717: call VecScale(kval,0.5d+0,ierr)
1719: call VecPointwiseDivide(steval,eval,rval,ierr)
1720: call VecWAXPY(seval,-1.0d+0,kval,steval,ierr)
1722: call VecCopy(seval,tval,ierr)
1723: call VecScale(tval,csubvi,ierr)
1725: c
1726: c plot tval vector
1727: c
1728: call VecView(tval,view3,ierr)
1729: c
1730: c make xmgr plot of tval
1731: c
1732: call FormXmgr(tval,4,ierr)
1734: c
1735: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1736: c
1738: call VecPointwiseMult(rval,seval,pval,ierr)
1739: call VecScale(pval,gm1,ierr)
1740: c
1741: c plot pval vector
1742: c
1743: call VecView(pval,view4,ierr)
1744: c
1745: c make xmgr plot of pval
1746: c
1747: call FormXmgr(pval,5,ierr)
1748: c
1749: c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1750: c
1756: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1757: ! Free work space. All PETSc objects should be destroyed when they
1758: ! are no longer needed.
1759: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1761: call VecDestroy(rval, ierr)
1762: call VecDestroy(uval, ierr)
1763: call VecDestroy(ruval,ierr)
1764: call VecDestroy(eval, ierr)
1765: call VecDestroy(seval, ierr)
1766: call VecDestroy(pval, ierr)
1767: call VecDestroy(kval, ierr)
1768: call VecDestroy(tval, ierr)
1769: call VecDestroy(steval, ierr)
1771: call ISDestroy(rfrom, ierr)
1772: call ISDestroy(rto, ierr)
1774: call ISDestroy(rufrom,ierr)
1775: call ISDestroy(ruto, ierr)
1777: call ISDestroy(efrom, ierr)
1778: call ISDestroy(eto, ierr)
1781: return
1782: end
1783: subroutine FormInitialGuess(x,ierr)
1784: ! ---------------------------------------------------------------------
1785: !
1786: ! FormInitialGuess - Forms initial approximation.
1787: !
1788: ! Input Parameter:
1789: ! x - vector
1790: !
1791: ! Output Parameters:
1792: ! x - vector
1793: ! ierr - error code
1794: !
1795: ! Notes:
1796: ! This routine serves as a wrapper for the lower-level routine
1797: ! "ApplicationInitialGuess", where the actual computations are
1798: ! done using the standard Fortran style of treating the local
1799: ! vector data as a multidimensional array over the local mesh.
1800: ! This routine merely accesses the local vector data via
1801: ! VecGetArray() and VecRestoreArray().
1802: !
1803: implicit none
1805: #include include/finclude/petsc.h
1806: #include include/finclude/petscvec.h
1807: #include include/finclude/petscmat.h
1808: #include include/finclude/petscsnes.h
1810: ! Input/output variables:
1811: Vec x
1812: integer ierr
1814: ! Declarations for use with local arrays:
1815: PetscScalar lx_v(0:1)
1816: PetscOffset lx_i
1818: 0
1820: c
1821: c get a pointer to x
1822: c
1823: call VecGetArray(x,lx_v,lx_i,ierr)
1824: c
1825: c Compute initial guess
1826: c
1827: call ApplicationInitialGuess(lx_v(lx_i),ierr)
1828: c
1829: c Restore vector x
1830: c
1831: call VecRestoreArray(x,lx_v,lx_i,ierr)
1833: return
1834: end
1835: subroutine FormXmgr(x,ivar,ierr)
1836: ! ---------------------------------------------------------------------
1837: !
1838: ! FormXmgr - Forms Xmgr output
1839: !