Actual source code: ex37f90.F90
1: #include "finclude/petscdef.h"
2: !
3: ! Notes:
4: ! This uses Fortran 90 free-form, this means the lines can be up to 132 columns wide
5: ! CHKR(ierr) is put on the same line as the call except when it violates the 132 columns
6: !
7: ! Eventually this file can be split into several files
8: !
9: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10: ! Error handler that aborts when error is detected
11: !
12: subroutine HE1(ierr,line)
13: use mex37f90
14: use mex37f90interfaces
16: call PetscError(ierr,line,0,'',ierr)
17: call MPI_Abort(PETSC_COMM_WORLD,ierr,ierr)
18: return
19: end
20: #define CHKQ(n) if(n .ne. 0)call HE1(n,__LINE__)
22: ! Error handler forces traceback of where error occurred
23: !
24: subroutine HE2(ierr,line)
25: use mex37f90
26: use mex37f90interfaces
28: call PetscError(ierr,line,0,'',ierr)
29: return
30: end
31: #define CHKR(n) if(n .ne. 0)then;call HE2(n,__LINE__);return;endif
33: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34: ! Utilities to map between global (linear algebra) to local (ghosted, physics) representations
35: !
36: ! Allocates the local work arrays and vectors
37: subroutine LocalFormCreate(dm,lf,ierr)
38: use mex37f90
39: use mex37f90interfaces
40: implicit none
41: DMComposite dm
42: type(LocalForm) lf
43: PetscErrorCode ierr
45: call DMCompositeGetEntries(dm,lf%da,lf%np,lf%da,lf%np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
46: call DAGetLocalInfof90(lf%da,lf%dainfo,ierr);CHKR(ierr)
47: call DMCompositeGetLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
48: CHKR(ierr)
49: return
50: end
52: ! Updates the (ghosted) IHX and Core arrays from the global vector x
53: subroutine LocalFormUpdate(dm,x,lf,ierr)
54: use mex37f90
55: use mex37f90interfaces
56: implicit none
57: DMComposite dm
58: type(LocalForm) lf
59: PetscErrorCode ierr
60: Vec x
62: call DMCompositeScatter(dm,x,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
63: CHKR(ierr)
64: call DAVecGetArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
65: call DAVecGetArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
66: return
67: end
69: ! You should call this when you no longer need access to %IHX and %Core
70: subroutine LocalFormRestore(dm,lf,ierr)
71: use mex37f90
72: use mex37f90interfaces
73: implicit none
74: DMComposite dm
75: type(LocalForm) lf
76: PetscErrorCode ierr
78: call DAVecRestoreArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
79: call DAVecRestoreArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
80: return
81: end
83: ! Destroys the data structure
84: subroutine LocalFormDestroy(dm,lf,ierr)
85: use mex37f90
86: use mex37f90interfaces
87: implicit none
88: DMComposite dm
89: type(LocalForm) lf
90: PetscErrorCode ierr
92: call DMCompositeRestoreLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
93: CHKR(ierr)
94: return
95: end
97: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
98: ! Implements a nonlinear solver for a simple domain
99: ! consisting of 2 zero dimensional problems linked by
100: ! 2 one dimensional problems.
101: !
102: ! Channel1
103: ! -------------------------
104: ! Pool 1 Pool 2
105: ! -------------------------
106: ! Channel2
107: !VAM
108: !VAM
109: !VAM
110: !VAM Hot Pool 1
111: !VAM | |
112: !VAM | |
113: !VAM | |
114: !VAM | |
115: !VAM Core Channel 2 | | IHX Channel 1
116: !VAM | |
117: !VAM | |
118: !VAM | |
119: !VAM | |
120: !VAM Cold Pool 2
121: !VAM
122: !
123: ! For Newton's method with block Jacobi (one block for
124: ! each channel and one block for each pool) run with
125: !
126: ! -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
127: ! -help lists all options
129: program ex37f90
130: use mex37f90
131: #include finclude/petscsys.h
132: #include finclude/petscviewer.h
133: #include finclude/petscvec.h
135: DMMGArray dmmg
136: DMMG dmmglevel
137: PetscErrorCode ierr
138: DA da
139: DMComposite dm
140: type(appctx) app
141: external FormInitialGuess,FormFunction,FormGraph
142: Vec x
144: double precision time
145: integer i
146: !BARRY
147: PetscViewer view0,view1
148: !BARRY
150: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
152: ! Create the composite DM object that manages the grid
154: call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
155: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
156: CHKR(ierr)
157: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
158: call DADestroy(da,ierr);CHKR(ierr)
159: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
160: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
161: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
163: ! Create solver object and associate it with the unknowns (on the grid)
165: call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
166: call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
167: call DMCompositeDestroy(dm,ierr);CHKR(ierr)
168: call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr) ! currently only one level solver
169: call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
170: call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(ierr)
171: call DMMGSetFromOptions(dmmg,ierr)
172: !BARRY
173: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'core',0,0,300,300,view0,ierr)
174: CHKR(ierr)
175: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'ihx',320,0,300,300,view1,ierr)
176: CHKR(ierr)
177: !BARRY
179: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
180: call VecDuplicate(x,app%xold,ierr);CHKR(ierr)
181: call VecCopy(x,app%xold,ierr);CHKR(ierr)
182: !BARRY
183: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
184: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
185: !BARRY
187: time = 0.0d+0
189: write(*,*)
190: write(*,*) 'initial time'
191: write(*,*)
193: do i = 1,app%nstep
195: time = time + app%dt
197: write(*,*)
198: write(*,*) 'time =',time
199: write(*,*)
200: !
201: ! copy new to old
202: !
203: !
204: ! Solve the nonlinear system
205: !
206: call DMMGSolve(dmmg,ierr);CHKR(ierr)
208: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
209: call VecCopy(x,app%xold,ierr);CHKR(ierr)
210: !BARRY
211: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
212: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
213: !BARRY
214: enddo
215: !BARRY
216: call PetscViewerDestroy(view0,ierr);CHKR(ierr)
217: call PetscViewerDestroy(view1,ierr);CHKR(ierr)
218: !BARRY
219: write(*,*)
220: write(*,*) 'final time'
221: write(*,*)
223: call VecDestroy(app%xold,ierr);CHKR(ierr)
224: call DMMGDestroy(dmmg,ierr);CHKR(ierr)
225: call PetscFinalize(ierr)
226: end
228: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229: !
230: ! The physics :-)
232: subroutine FormFunction(snes,x,f,dmmg,ierr)
233: use mex37f90
234: use mex37f90interfaces
236: SNES snes
237: Vec x,f
238: DMMG dmmg
239: PetscErrorCode ierr
241: DMComposite dm
242: Vec fvc1,fvc2
243: PetscInt np,i
244: DA da
245: type(poolfield), pointer :: fHotPool,fColdPool
246: type(channelfield), pointer :: fIHX(:),fCore(:)
247: type(appctx), pointer :: app
248: PetscMPIInt rank
249: type(LocalForm) xlf,xoldlf
251: double precision dhhpl, mult, dhcpl, phpco, pcpio, pcpci, phpii
253: logical debug
255: debug = .false.
257: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
259: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr) ! access the grid information
261: call LocalFormCreate(dm,xlf,ierr) ! Access the local parts of x
262: call LocalFormUpdate(dm,x,xlf,ierr)
264: ! Access the global (non-ghosted) parts of f
265: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
266: call DMCompositeGetAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
267: call DAVecGetArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
268: call DAVecGetArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
270: !BARRY
271: !
272: ! At this point I need to create old time values from "xold" passed in through
273: ! app for
274: !
275: ! xHotPool%vol, xHotPool%vel, xColdPool%vol, xColdPool%vel
276: ! xIHX(i)%press, xCore(i)%press
277: !
278: ! I don't know how to do this.
279: !
280: !BARRY
281: call LocalFormCreate(dm,xoldlf,ierr) ! Access the local parts of xold
282: call LocalFormUpdate(dm,app%xold,xoldlf,ierr)
284: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
285: ! fPool only lives on zeroth process, ghosted xPool lives on all processes
286: if (rank == 0) then
287: !
288: ! compute velocity into hotpool from core
289: !
290: dhhpl = app%dhhpl0 + ( (xlf%HotPool%vol - app%hpvol0) / app%ahp )
291: phpco = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhco) )
292: mult = app%dt / (app%dxc * app%rho)
293: fHotPool%vel = xlf%HotPool%vel - (mult * (xlf%Core(app%nxc-1)%press - phpco) ) + (abs(xlf%HotPool%vel) * xlf%HotPool%vel)
294: !
295: ! compute change in hot pool volume
296: !
297: fHotPool%vol = xlf%HotPool%vol - app%hpvol0 - (app%dt * app%acore * (xlf%HotPool%vel -xlf%ColdPool%vel) )
298: !
299: ! compute velocity into coldpool from IHX
300: !
301: dhcpl = app%dhcpl0 + ( (xlf%ColdPool%vol - app%cpvol0) / app%acp )
302: pcpio = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhio) )
303: mult = app%dt / (app%dxc * app%rho)
304: fColdPool%vel = xlf%ColdPool%vel-(mult*(xlf%IHX(app%nxc-1)%press-pcpio))+(abs(xlf%ColdPool%vel)*xlf%ColdPool%vel)
305: !
306: ! compute the change in cold pool volume
307: !
308: fColdPool%vol = xlf%ColdPool%vol - app%cpvol0 - (app%dt * app%aihx * (xlf%ColdPool%vel - xlf%HotPool%vel) )
309: endif
310: !
311: ! Compute the pressure distribution in IHX and core
312: !
313: pcpci = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhci) )
314: phpii = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhii) )
316: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
318: fIHX(i)%press = xlf%IHX(i)%press - phpii - (app%rho * app%grav * dble(i) * app%dxi)
319: fCore(i)%press = xlf%Core(i)%press - pcpci + (app%rho * app%grav * dble(i) * app%dxc)
321: enddo
323: if (debug) then
324: write(*,*)
325: write(*,*) 'Hot vel,vol ',xlf%HotPool%vel,xlf%HotPool%vol
326: write(*,*) 'delta p = ', xlf%Core(app%nxc-1)%press - phpco,xlf%Core(app%nxc-1)%press,phpco
327: write(*,*)
329: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
330: write(*,*) 'xlf%IHX(',i,')%press = ',xlf%IHX(i)%press
331: enddo
333: write(*,*)
334: write(*,*) 'Cold vel,vol ',xlf%ColdPool%vel,xlf%ColdPool%vol
335: write(*,*) 'delta p = ',xlf%IHX(app%nxc-1)%press - pcpio,xlf%IHX(app%nxc-1)%press, pcpio
336: write(*,*)
339: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
340: write(*,*) 'xlf%Core(',i,')%press = ',xlf%Core(i)%press
341: enddo
343: endif
345: call DAVecRestoreArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
346: call DAVecRestoreArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
347: call DMCompositeRestoreAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
349: call LocalFormRestore(dm,xoldlf,ierr)
350: call LocalFormDestroy(dm,xoldlf,ierr)
352: call LocalFormRestore(dm,xlf,ierr)
353: call LocalFormDestroy(dm,xlf,ierr)
354: return
355: end
357: subroutine FormGraph(dmmg,x,view0,view1,ierr)
359: ! --------------------------------------------------------------------------------------
360: !
361: ! FormGraph - Forms Graphics output
362: !
363: ! Input Parameter:
364: ! x - vector
365: ! time - scalar
366: !
367: ! Output Parameters:
368: ! ierr - error code
369: !
370: ! Notes:
371: ! This routine serves as a wrapper for the lower-level routine
372: ! "ApplicationXmgr", where the actual computations are
373: ! done using the standard Fortran style of treating the local
374: ! vector data as a multidimensional array over the local mesh.
375: ! This routine merely accesses the local vector data via
376: ! VecGetArray() and VecRestoreArray().
377: !
378: use mex37f90
379: use mex37f90interfaces
381: Vec x,xvc1,xvc2 !,corep,ihxp
382: DMMG dmmg
383: PetscErrorCode ierr
384: DMComposite dm
385: PetscInt np !,i
386: DA da
387: type(DALocalInfof90) dainfo
388: type(poolfield), pointer :: xHotPool,xColdPool
389: type(channelfield), pointer :: xIHX(:),xCore(:)
390: type(appctx), pointer :: app
392: PetscViewer view0,view1
394: integer iplotnum
395: save iplotnum
396: character*8 grfile
398: data iplotnum / -1 /
399: !BARRY
400: !
401: ! This is my attempt to get the information out of vector "x" to plot
402: ! it. I need to be able to build xCore(i)%press and xIHX(i)%press
403: ! from the vector "x" and I do not know how to do this
404: !
405: !BARRY
407: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
408: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr) ! Access the grid information
410: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
411: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
413: call DMCompositeGetLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
414: call DMCompositeScatter(dm,x,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
415: call DAVecGetArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
416: call DAVecGetArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
418: !
419: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
420: !
421: iplotnum = iplotnum + 1
422: 0
423: !
424: !
425: ! plot corep vector
426: !
427: call VecView(xvc1,view0,ierr);CHKR(ierr)
428: !
429: ! make xmgr plot of corep
430: !
431: write(grfile,4000) iplotnum
432: 4000 format('CoreP',i3.3)
434: open (unit=44,file=grfile,status='unknown')
436: do i = 1,app%nxc
437: write(44,1000) dble(i)*app%dxc, xCore(i)%press
438: enddo
440: close(44)
441: !
442: ! plot ihxp vector
443: !
444: call VecView(xvc2,view1,ierr);CHKR(ierr)
445: !
446: ! make xmgr plot of ihxp
447: !
448: write(grfile,8000) iplotnum
449: 8000 format('IHXPr',i3.3)
451: open (unit=44,file=grfile,status='unknown')
453: do i = 1,app%nxc
454: write(44,1000) dble(i)*app%dxi, xIHX(i)%press
455: enddo
457: close(44)
461: 1000 format(3(e18.12,2x))
463: call DAVecRestoreArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
464: call DAVecRestoreArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
465: call DMCompositeRestoreLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
467: return
468: end
470: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
472: subroutine FormInitialGuess(dmmg,v,ierr)
473: use mex37f90
474: use mex37f90interfaces
476: DMMG dmmg
477: Vec v
478: PetscErrorCode ierr
480: DMComposite dm
481: Vec vc1,vc2
482: PetscInt np,i
483: DA da
484: type(DALocalInfof90) dainfo
485: type(poolfield), pointer :: HotPool,ColdPool
486: type(channelfield), pointer :: IHX(:),Core(:)
487: type(appctx), pointer :: app
488: PetscMPIInt rank
490: double precision pcpci, phpii
492: logical debug
494: debug = .false.
496: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
498: ! Access the grid information
500: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
501: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
502: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
504: ! Acess the vector unknowns in global (nonghosted) form
506: call DMCompositeGetAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
507: call DAVecGetArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
508: call DAVecGetArrayf90(da,vc2,Core,ierr);CHKR(ierr)
510: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
511: !
512: ! Set the input values
513: !
515: app%P0 = 1.0d+5
516: app%rho = 1.0d+3
517: app%grav = 9.8d+0
519: app%dhhpl0 = 12.2d+0
520: app%dhcpl0 = 10.16d+0
522: app%lenc = 3.0d+0
523: app%leni = 4.46d+0
525: app%dhci = 0.83d+0
526: app%dhco = app%dhci + app%lenc
528: app%dhii = 7.83d+0
529: app%dhio = app%dhii - app%leni
531: app%dxc = app%lenc / dble(app%nxc)
532: app%dxi = app%leni / dble(app%nxc)
534: app%dt = 1.0d+0
536: app%ahp = 7.0d+0
537: app%acp = 7.0d+0
539: app%acore = 0.8d+0
540: app%aihx = 5.0d-2
542: app%hpvol0 = app%dhhpl0 * app%ahp
543: app%cpvol0 = app%dhcpl0 * app%acp
544: !
545: ! the pools are only unknowns on process 0
546: !
547: if (rank == 0) then
548: HotPool%vel = -1.0d+0
549: HotPool%vol = app%hpvol0
550: ColdPool%vel = 1.0d+0
551: ColdPool%vol = app%cpvol0
552: endif
553: !
554: ! Construct and initial guess at the pressure
555: !
556: pcpci = app%P0 + ( app%rho * app%grav * (app%dhcpl0 - app%dhci) )
557: phpii = app%P0 + ( app%rho * app%grav * (app%dhhpl0 - app%dhii) )
559: if (debug) then
560: write(*,*)
561: write(*,*) 'pcpci = ',pcpci
562: write(*,*) 'phpii = ',phpii
563: write(*,*) 'app%P0 = ',app%P0
564: write(*,*) 'dhcpl0 - app%dhci ',app%dhcpl0 - app%dhci,app%dhcpl0, app%dhci
565: write(*,*) 'dhhpl0 - app%dhii ',app%dhhpl0 - app%dhii,app%dhhpl0, app%dhii
566: write(*,*)
567: endif
569: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
571: IHX(i)%press = phpii + (app%rho * app%grav * dble(i) * app%dxi)
573: Core(i)%press = pcpci - (app%rho * app%grav * dble(i) * app%dxc)
575: enddo
577: call DAVecRestoreArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
578: call DAVecRestoreArrayf90(da,vc2,Core,ierr);CHKR(ierr)
579: call DMCompositeRestoreAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
581: CHKMEMQ
582: return
583: end