Actual source code: ex35f90.F90
1: #include "finclude/petscdef.h"
3: ! Error handler that aborts when error is detected
4: !
5: subroutine HE1(ierr,line)
6: use mex35f90
7: use mex35f90interfaces
9: call PetscError(ierr,line,0,'',ierr)
10: call MPI_Abort(PETSC_COMM_WORLD,ierr,ierr)
11: return
12: end
13: #define CHKQ(n) if(n .ne. 0)call HE1(n,__LINE__)
15: ! Error handler forces traceback of where error occurred
16: !
17: subroutine HE2(ierr,line)
18: use mex35f90
19: use mex35f90interfaces
21: call PetscError(ierr,line,0,'',ierr)
22: return
23: end
24: #define CHKR(n) if(n .ne. 0)then;call HE2(n,__LINE__);return;endif
26: !
27: ! Implements a nonlinear solver for a simple domain
28: ! consisting of 2 zero dimensional problems linked by
29: ! 2 one dimensional problems.
30: !
31: ! Channel1
32: ! -------------------------
33: ! Pool 1 Pool 2
34: ! -------------------------
35: ! Channel2
36: !VAM
37: !VAM
38: !VAM
39: !VAM Hot Pool 1
40: !VAM | |
41: !VAM | |
42: !VAM | |
43: !VAM | |
44: !VAM Core Channel 2 | | IHX Channel 1
45: !VAM | |
46: !VAM | |
47: !VAM | |
48: !VAM | |
49: !VAM Cold Pool 2
50: !VAM
51: !
52: ! For Newton's method with block Jacobi (one block for
53: ! each channel and one block for each pool) run with
54: !
55: ! -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
56: ! -help lists all options
58: program ex35f90
59: use mex35f90
60: ! use mex35f90interfaces
61: #include finclude/petscsys.h
62: #include finclude/petscviewer.h
63: #include finclude/petscvec.h
66: DMMGArray dmmg
67: DMMG dmmglevel
68: PetscErrorCode ierr
69: DA da
70: DMComposite dm
71: type(appctx) app
72: external FormInitialGuess,FormFunction,FormGraph
73: Vec x
75: double precision time
76: integer i
77: !BARRY
78: PetscViewer view0,view1
79: !BARRY
81: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
83: ! Create the composite DM object that manages the grid
85: call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
86: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
87: CHKR(ierr)
88: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
89: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
90: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
91: call DADestroy(da,ierr);CHKR(ierr)
92: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
94: ! Create solver object and associate it with the unknowns (on the grid)
96: call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
97: call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
98: call DMCompositeDestroy(dm,ierr);CHKR(ierr)
99: call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr) ! currently only one level solver
100: call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
101: call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(ierr)
102: call DMMGSetFromOptions(dmmg,ierr)
103: !BARRY
104: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'core',0,0,300,300,view0,ierr)
105: CHKR(ierr)
106: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'ihx',320,0,300,300,view1,ierr)
107: CHKR(ierr)
108: !BARRY
110: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
111: call VecDuplicate(x,app%xold,ierr);CHKR(ierr)
112: call VecCopy(x,app%xold,ierr);CHKR(ierr)
113: !BARRY
114: write(*,*) 'Before call to FormGraph'
115: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
116: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
117: !BARRY
119: time = 0.0d+0
121: write(*,*)
122: write(*,*) 'initial time'
123: write(*,*)
125: do i = 1,app%nstep
127: time = time + app%dt
129: write(*,*)
130: write(*,*) 'time =',time
131: write(*,*)
132: !
133: ! copy new to old
134: !
135: !
136: ! Solve the nonlinear system
137: !
138: call DMMGSolve(dmmg,ierr);CHKR(ierr)
139: call DMMGSetInitialGuess(dmmg,PETSC_NULL_FUNCTION,ierr);CHKR(ierr)
141: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
142: call VecCopy(x,app%xold,ierr);CHKR(ierr)
143: !BARRY
144: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
145: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
146: !BARRY
147: enddo
148: !BARRY
149: call PetscViewerDestroy(view0,ierr);CHKR(ierr)
150: call PetscViewerDestroy(view1,ierr);CHKR(ierr)
151: !BARRY
152: write(*,*)
153: write(*,*) 'final time'
154: write(*,*)
156: call VecDestroy(app%xold,ierr);CHKR(ierr)
157: call DMMGDestroy(dmmg,ierr);CHKR(ierr)
158: call PetscFinalize(ierr)
159: end
161: subroutine FormFunction(snes,x,f,dmmg,ierr)
162: use mex35f90
163: use mex35f90interfaces
165: SNES snes
166: Vec x,f
167: DMMG dmmg
168: PetscErrorCode ierr
170: DMComposite dm
171: Vec fvc1,fvc2,xvc1,xvc2
172: PetscInt np,i
173: DA da
174: type(DALocalInfof90) dainfo
175: type(poolfield), pointer :: fHotPool,fColdPool
176: type(poolfield), pointer :: xHotPool,xColdPool
177: type(channelfield), pointer :: fIHX(:),fCore(:),xIHX(:),xCore(:)
178: type(appctx), pointer :: app
179: PetscMPIInt rank
181: double precision dhhpl, mult, dhcpl, phpco, pcpio, pcpci, phpii
183: logical debug
185: debug = .false.
187: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
189: ! Access the grid information
191: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
192: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
193: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
195: ! Access the local (ghosted) parts of x
197: call DMCompositeGetLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
198: call DMCompositeScatter(dm,x,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
200: call DAVecGetArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
201: call DAVecGetArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
203: ! Access the global (non-ghosted) parts of f
205: call DMCompositeGetAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
206: call DAVecGetArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
207: call DAVecGetArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
209: !BARRY
210: !
211: ! At this point I need to create old time values from "xold" passed in through
212: ! app for
213: !
214: ! xHotPool%vol, xHotPool%vel, xColdPool%vol, xColdPool%vel
215: ! xIHX(i)%press, xCore(i)%press
216: !
217: ! I don't know how to do this.
218: !
219: !BARRY
221: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
222: ! fPool only lives on zeroth process, ghosted xPool lives on all processes
223: if (rank == 0) then
224: !
225: ! compute velocity into hotpool from core
226: !
227: dhhpl = app%dhhpl0 + ( (xHotPool%vol - app%hpvol0) / app%ahp )
228: phpco = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhco) )
229: mult = app%dt / (app%dxc * app%rho)
230: fHotPool%vel = xHotPool%vel - (mult * (xCore(app%nxc-1)%press - phpco) ) + (abs(xHotPool%vel) * xHotPool%vel)
231: !
232: ! compute change in hot pool volume
233: !
234: fHotPool%vol = xHotPool%vol - app%hpvol0 - (app%dt * app%acore * (xHotPool%vel -xColdPool%vel) )
235: !
236: ! compute velocity into coldpool from IHX
237: !
238: dhcpl = app%dhcpl0 + ( (xColdPool%vol - app%cpvol0) / app%acp )
239: pcpio = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhio) )
240: mult = app%dt / (app%dxc * app%rho)
241: fColdPool%vel = xColdPool%vel - (mult * (xIHX(app%nxc-1)%press - pcpio) ) + (abs(xColdPool%vel) * xColdPool%vel)
242: !
243: ! compute the change in cold pool volume
244: !
245: fColdPool%vol = xColdPool%vol - app%cpvol0 - (app%dt * app%aihx * (xColdPool%vel - xHotPool%vel) )
246: endif
247: !
248: ! Compute the pressure distribution in IHX and core
249: !
250: pcpci = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhci) )
251: phpii = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhii) )
253: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
255: fIHX(i)%press = xIHX(i)%press - phpii - (app%rho * app%grav * dble(i) * app%dxi)
257: fCore(i)%press = xCore(i)%press - pcpci + (app%rho * app%grav * dble(i) * app%dxc)
259: enddo
261: if (debug) then
262: write(*,*)
263: write(*,*) 'Hot vel,vol ',xHotPool%vel,xHotPool%vol
264: write(*,*) 'delta p = ', xCore(app%nxc-1)%press - phpco,xCore(app%nxc-1)%press,phpco
265: write(*,*)
267: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
268: write(*,*) 'xIHX(',i,')%press = ',xIHX(i)%press
269: enddo
271: write(*,*)
272: write(*,*) 'Cold vel,vol ',xColdPool%vel,xColdPool%vol
273: write(*,*) 'delta p = ',xIHX(app%nxc-1)%press - pcpio,xIHX(app%nxc-1)%press, pcpio
274: write(*,*)
277: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
278: write(*,*) 'xCore(',i,')%press = ',xCore(i)%press
279: enddo
281: endif
283: call DAVecRestoreArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
284: call DAVecRestoreArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
285: call DMCompositeRestoreLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
287: call DAVecRestoreArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
288: call DAVecRestoreArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
289: call DMCompositeRestoreAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
291: return
292: end
293: subroutine FormGraph(dmmg,x,view0,view1,ierr)
294: ! ---------------------------------------------------------------------
295: !
296: ! FormGraph - Forms Graphics output
297: !
298: ! Input Parameter:
299: ! x - vector
300: ! time - scalar
301: !
302: ! Output Parameters:
303: ! ierr - error code
304: !
305: ! Notes:
306: ! This routine serves as a wrapper for the lower-level routine
307: ! "ApplicationXmgr", where the actual computations are
308: ! done using the standard Fortran style of treating the local
309: ! vector data as a multidimensional array over the local mesh.
310: ! This routine merely accesses the local vector data via
311: ! VecGetArray() and VecRestoreArray().
312: !
313: use mex35f90
314: use mex35f90interfaces
316: Vec x,xvc1,xvc2 !,corep,ihxp
317: DMMG dmmg
318: PetscErrorCode ierr
319: DMComposite dm
320: PetscInt np !,i
321: DA da
322: type(DALocalInfof90) dainfo
323: type(poolfield), pointer :: HotPool,ColdPool
324: type(poolfield), pointer :: xHotPool,xColdPool
325: type(channelfield), pointer :: xIHX(:),xCore(:)
326: type(appctx), pointer :: app
327: PetscMPIInt rank
329: PetscViewer view0,view1
331: integer iplotnum
332: save iplotnum
333: character*8 grfile
335: data iplotnum / -1 /
336: !BARRY
337: !
338: ! This is my attempt to get the information out of vector "x" to plot
339: ! it. I need to be able to build xCore(i)%press and xIHX(i)%press
340: ! from the vector "x" and I do not know how to do this
341: !
342: !BARRY
344: write(*,*)
345: write(*,*) 'inside of FormGraph'
346: write(*,*)
348: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
350: write(*,*)
351: write(*,*) 'after DMMGGetUser'
352: write(*,*)
354: ! Access the grid information
356: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
358: write(*,*)
359: write(*,*) 'after DMMGGetDM'
360: write(*,*)
362: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
363: write(*,*)
364: write(*,*) 'after DMCompositeGetEntries'
365: write(*,*)
366: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
367: write(*,*)
368: write(*,*) 'after DAGetLocalInfof90'
369: write(*,*)
370: !BARRY
371: !
372: ! I think that the code dies in the call below.
373: !
374: !BARRY
375: call DMCompositeGetLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
376: call DMCompositeScatter(dm,x,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
377: call DAVecGetArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
378: write(*,*)
379: write(*,*) 'after DAVecGetArrayf90(da,xvc1,xIHX,ierr)'
380: write(*,*)
381: call DAVecGetArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
382: write(*,*)
383: write(*,*) 'after DAVecGetArrayf90(da,xvc2,xCore,ierr)'
384: write(*,*)
387: !
388: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
389: !
390: iplotnum = iplotnum + 1
391: 0
392: !
393: !
394: ! plot corep vector
395: !
396: call VecView(xvc1,view0,ierr);CHKR(ierr)
397: !
398: ! make xmgr plot of corep
399: !
400: write(grfile,4000) iplotnum
401: 4000 format('CoreP',i3.3)
403: open (unit=44,file=grfile,status='unknown')
405: do i = 1,app%nxc
406: write(44,1000) dble(i)*app%dxc, xCore(i)%press
407: enddo
409: close(44)
410: !
411: ! plot ihxp vector
412: !
413: call VecView(xvc2,view1,ierr);CHKR(ierr)
414: !
415: ! make xmgr plot of ihxp
416: !
417: write(grfile,8000) iplotnum
418: 8000 format('IHXPr',i3.3)
420: open (unit=44,file=grfile,status='unknown')
422: do i = 1,app%nxc
423: write(44,1000) dble(i)*app%dxi, xIHX(i)%press
424: enddo
426: close(44)
430: 1000 format(3(e18.12,2x))
432: call DAVecRestoreArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
433: call DAVecRestoreArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
434: call DMCompositeRestoreLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
435: return
436: end
439: subroutine FormInitialGuess(dmmg,v,ierr)
440: use mex35f90
441: use mex35f90interfaces
443: DMMG dmmg
444: Vec v
445: PetscErrorCode ierr
447: DMComposite dm
448: Vec vc1,vc2
449: PetscInt np,i
450: DA da
451: type(DALocalInfof90) dainfo
452: type(poolfield), pointer :: HotPool,ColdPool
453: type(channelfield), pointer :: IHX(:),Core(:)
454: type(appctx), pointer :: app
455: PetscMPIInt rank
457: double precision pcpci, phpii
459: logical debug
461: debug = .false.
463: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
465: ! Access the grid information
467: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
468: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
469: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
471: ! Acess the vector unknowns in global (nonghosted) form
473: call DMCompositeGetAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
474: call DAVecGetArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
475: call DAVecGetArrayf90(da,vc2,Core,ierr);CHKR(ierr)
477: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
478: !
479: ! Set the input values
480: !
482: app%P0 = 1.0d+5
483: app%rho = 1.0d+3
484: app%grav = 9.8d+0
486: app%dhhpl0 = 12.2d+0
487: app%dhcpl0 = 10.16d+0
489: app%lenc = 3.0d+0
490: app%leni = 4.46d+0
492: app%dhci = 0.83d+0
493: app%dhco = app%dhci + app%lenc
495: app%dhii = 7.83d+0
496: app%dhio = app%dhii - app%leni
498: app%dxc = app%lenc / dble(app%nxc)
499: app%dxi = app%leni / dble(app%nxc)
501: app%dt = 1.0d+0
503: app%ahp = 7.0d+0
504: app%acp = 7.0d+0
506: app%acore = 0.8d+0
507: app%aihx = 5.0d-2
509: app%hpvol0 = app%dhhpl0 * app%ahp
510: app%cpvol0 = app%dhcpl0 * app%acp
511: !
512: ! the pools are only unknowns on process 0
513: !
514: if (rank == 0) then
515: HotPool%vel = -1.0d+0
516: HotPool%vol = app%hpvol0
517: ColdPool%vel = 1.0d+0
518: ColdPool%vol = app%cpvol0
519: endif
520: !
521: ! Construct and initial guess at the pressure
522: !
523: pcpci = app%P0 + ( app%rho * app%grav * (app%dhcpl0 - app%dhci) )
524: phpii = app%P0 + ( app%rho * app%grav * (app%dhhpl0 - app%dhii) )
526: if (debug) then
527: write(*,*)
528: write(*,*) 'pcpci = ',pcpci
529: write(*,*) 'phpii = ',phpii
530: write(*,*) 'app%P0 = ',app%P0
531: write(*,*) 'dhcpl0 - app%dhci ',app%dhcpl0 - app%dhci,app%dhcpl0, app%dhci
532: write(*,*) 'dhhpl0 - app%dhii ',app%dhhpl0 - app%dhii,app%dhhpl0, app%dhii
533: write(*,*)
534: endif
536: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
538: IHX(i)%press = phpii + (app%rho * app%grav * dble(i) * app%dxi)
540: Core(i)%press = pcpci - (app%rho * app%grav * dble(i) * app%dxc)
542: enddo
544: call DAVecRestoreArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
545: call DAVecRestoreArrayf90(da,vc2,Core,ierr);CHKR(ierr)
546: call DMCompositeRestoreAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
548: CHKMEMQ
549: return
550: end