Actual source code: ex37f90.F90
1: #define PETSC_AVOID_DECLARATIONS
2: #include include/finclude/petscall.h
3: !
4: ! Notes:
5: ! This uses Fortran 90 free-form, this means the lines can be up to 132 columns wide
6: ! CHKR(ierr) is put on the same line as the call except when it violates the 132 columns
7: !
8: ! Eventually this file can be split into several files
9: !
10: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11: ! Error handler that aborts when error is detected
12: !
13: subroutine HE1(ierr,line)
14: use mex37f90
15: use mex37f90interfaces
17: call PetscError(ierr,line,0,'',ierr)
18: call MPI_Abort(PETSC_COMM_WORLD,ierr,ierr)
19: return
20: end
21: #define CHKQ(n) if(n .ne. 0)call HE1(n,__LINE__)
23: ! Error handler forces traceback of where error occurred
24: !
25: subroutine HE2(ierr,line)
26: use mex37f90
27: use mex37f90interfaces
29: call PetscError(ierr,line,0,'',ierr)
30: return
31: end
32: #define CHKR(n) if(n .ne. 0)then;call HE2(n,__LINE__);return;endif
34: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35: ! Utilities to map between global (linear algebra) to local (ghosted, physics) representations
36: !
37: ! Allocates the local work arrays and vectors
38: subroutine LocalFormCreate(dm,lf,ierr)
39: use mex37f90
40: use mex37f90interfaces
41: implicit none
42: DMComposite dm
43: type(LocalForm) lf
44: PetscErrorCode ierr
46: call DMCompositeGetEntries(dm,lf%da,lf%np,lf%da,lf%np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
47: call DAGetLocalInfof90(lf%da,lf%dainfo,ierr);CHKR(ierr)
48: call DMCompositeGetLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
49: CHKR(ierr)
50: return
51: end
53: ! Updates the (ghosted) IHX and Core arrays from the global vector x
54: subroutine LocalFormUpdate(dm,x,lf,ierr)
55: use mex37f90
56: use mex37f90interfaces
57: implicit none
58: DMComposite dm
59: type(LocalForm) lf
60: PetscErrorCode ierr
61: Vec x
63: call DMCompositeScatter(dm,x,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
64: CHKR(ierr)
65: call DAVecGetArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
66: call DAVecGetArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
67: return
68: end
70: ! You should call this when you no longer need access to %IHX and %Core
71: subroutine LocalFormRestore(dm,lf,ierr)
72: use mex37f90
73: use mex37f90interfaces
74: implicit none
75: DMComposite dm
76: type(LocalForm) lf
77: PetscErrorCode ierr
79: call DAVecRestoreArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
80: call DAVecRestoreArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
81: return
82: end
84: ! Destroys the data structure
85: subroutine LocalFormDestroy(dm,lf,ierr)
86: use mex37f90
87: use mex37f90interfaces
88: implicit none
89: DMComposite dm
90: type(LocalForm) lf
91: PetscErrorCode ierr
93: call DMCompositeRestoreLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
94: CHKR(ierr)
95: return
96: end
98: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99: ! Implements a nonlinear solver for a simple domain
100: ! consisting of 2 zero dimensional problems linked by
101: ! 2 one dimensional problems.
102: !
103: ! Channel1
104: ! -------------------------
105: ! Pool 1 Pool 2
106: ! -------------------------
107: ! Channel2
108: !VAM
109: !VAM
110: !VAM
111: !VAM Hot Pool 1
112: !VAM | |
113: !VAM | |
114: !VAM | |
115: !VAM | |
116: !VAM Core Channel 2 | | IHX Channel 1
117: !VAM | |
118: !VAM | |
119: !VAM | |
120: !VAM | |
121: !VAM Cold Pool 2
122: !VAM
123: !
124: ! For Newton's method with block Jacobi (one block for
125: ! each channel and one block for each pool) run with
126: !
127: ! -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
128: ! -help lists all options
130: program ex37f90
131: use mex37f90
132: #include include/finclude/petsc.h
133: #include include/finclude/petscviewer.h
134: #include include/finclude/petscvec.h
136: DMMGArray dmmg
137: DMMG dmmglevel
138: PetscErrorCode ierr
139: DA da
140: DMComposite dm
141: type(appctx) app
142: external FormInitialGuess,FormFunction,FormGraph
143: Vec x
145: double precision time
146: integer i
147: !BARRY
148: PetscViewer view0,view1
149: !BARRY
151: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
153: ! Create the composite DM object that manages the grid
155: call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
156: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
157: CHKR(ierr)
158: call DMCompositeAddDA(dm,da,ierr);CHKR(ierr)
159: call DADestroy(da,ierr);CHKR(ierr)
160: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
161: call DMCompositeAddDA(dm,da,ierr);CHKR(ierr)
162: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
164: ! Create solver object and associate it with the unknowns (on the grid)
166: call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
167: call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
168: call DMCompositeDestroy(dm,ierr);CHKR(ierr)
169: call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr) ! currently only one level solver
170: call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
171: call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(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