Actual source code: ex36f90.F90
1: #include "finclude/petscalldef.h"
3: ! Error handler that aborts when error is detected
4: !
5: subroutine HE1(ierr,line)
6: use mex36f90
7: use mex36f90interfaces
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 mex36f90
19: use mex36f90interfaces
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: ! Utilities to map between global (linear algebra) to local (ghosted, physics) representations
28: !
29: ! Allocates the local work arrays and vectors
30: subroutine LocalFormCreate(dm,lf,ierr)
31: use mex36f90
32: use mex36f90interfaces
33: implicit none
34: DMComposite dm
35: type(LocalForm) lf
36: PetscErrorCode ierr
38: call DMCompositeGetEntries(dm,lf%da,lf%np,lf%da,lf%np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
39: call DAGetLocalInfof90(lf%da,lf%dainfo,ierr);CHKR(ierr)
41: call DMCompositeGetLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
42: CHKR(ierr)
43: return
44: end
46: ! Updates the (ghosted) IHX and Core arrays from the global vector x
47: subroutine LocalFormUpdate(dm,x,lf,ierr)
48: use mex36f90
49: use mex36f90interfaces
50: implicit none
51: DMComposite dm
52: type(LocalForm) lf
53: PetscErrorCode ierr
54: Vec x
56: call DMCompositeScatter(dm,x,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
57: CHKR(ierr)
58: call DAVecGetArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
59: call DAVecGetArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
60: return
61: end
63: ! You should call this when you no longer need access to %IHX and %Core
64: subroutine LocalFormRestore(dm,lf,ierr)
65: use mex36f90
66: use mex36f90interfaces
67: implicit none
68: DMComposite dm
69: type(LocalForm) lf
70: PetscErrorCode ierr
72: call DAVecRestoreArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
73: call DAVecRestoreArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
74: return
75: end
77: ! Destroys the data structure
78: subroutine LocalFormDestroy(dm,lf,ierr)
79: use mex36f90
80: use mex36f90interfaces
81: implicit none
82: DMComposite dm
83: type(LocalForm) lf
84: PetscErrorCode ierr
86: call DMCompositeRestoreLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
87: CHKR(ierr)
88: return
89: end
91: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
92: ! Implements a nonlinear solver for a simple domain
93: ! consisting of 2 zero dimensional problems linked by
94: ! 2 one dimensional problems.
95: !
96: ! Channel1
97: ! -------------------------
98: ! Pool 1 Pool 2
99: ! -------------------------
100: ! Channel2
101: !VAM
102: !VAM
103: !VAM
104: !VAM Hot Pool 1
105: !VAM | |
106: !VAM | |
107: !VAM | |
108: !VAM | |
109: !VAM Core Channel 2 | | IHX Channel 1
110: !VAM | |
111: !VAM | |
112: !VAM | |
113: !VAM | |
114: !VAM Cold Pool 2
115: !VAM
116: !
117: ! For Newton's method with block Jacobi (one block for
118: ! each channel and one block for each pool) run with
119: !
120: ! -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
121: ! -help lists all options
123: program ex36f90
124: use mex36f90
125: ! use mex36f90interfaces
126: #include finclude/petsc.h
127: #include finclude/petscviewer.h
128: #include finclude/petscvec.h
131: DMMGArray dmmg
132: DMMG dmmglevel
133: PetscErrorCode ierr
134: DA da
135: DMComposite dm
136: type(appctx) app
137: external FormInitialGuess,FormFunction,FormGraph
138: Vec x
140: double precision time
141: integer i
142: !BARRY
143: PetscViewer view0,view1
144: !BARRY
146: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
148: ! Create the composite DM object that manages the grid
150: call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
151: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
152: CHKR(ierr)
153: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
154: call DADestroy(da,ierr);CHKR(ierr)
155: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
156: call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
157: call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
159: ! Create solver object and associate it with the unknowns (on the grid)
161: call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
162: call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
163: call DMCompositeDestroy(dm,ierr);CHKR(ierr)
164: call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr) ! currently only one level solver
165: call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
166: call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(ierr)
167: call DMMGSetFromOptions(dmmg,ierr)
168: !BARRY
169: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'core',0,0,300,300,view0,ierr)
170: CHKR(ierr)
171: call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'ihx',320,0,300,300,view1,ierr)
172: CHKR(ierr)
173: !BARRY
175: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
176: call VecDuplicate(x,app%xold,ierr);CHKR(ierr)
177: call VecCopy(x,app%xold,ierr);CHKR(ierr)
178: !BARRY
179: write(*,*) 'Before call to FormGraph'
180: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
181: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
182: !BARRY
184: time = 0.0d+0
186: write(*,*)
187: write(*,*) 'initial time'
188: write(*,*)
190: do i = 1,app%nstep
192: time = time + app%dt
194: write(*,*)
195: write(*,*) 'time =',time
196: write(*,*)
197: !
198: ! copy new to old
199: !
200: !
201: ! Solve the nonlinear system
202: !
203: call DMMGSolve(dmmg,ierr);CHKR(ierr)
205: call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
206: call VecCopy(x,app%xold,ierr);CHKR(ierr)
207: !BARRY
208: call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
209: call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
210: !BARRY
211: enddo
212: !BARRY
213: call PetscViewerDestroy(view0,ierr);CHKR(ierr)
214: call PetscViewerDestroy(view1,ierr);CHKR(ierr)
215: !BARRY
216: write(*,*)
217: write(*,*) 'final time'
218: write(*,*)
220: call VecDestroy(app%xold,ierr);CHKR(ierr)
221: call DMMGDestroy(dmmg,ierr);CHKR(ierr)
222: call PetscFinalize(ierr)
223: end
225: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226: !
227: ! The physics
229: subroutine FormFunction(snes,x,f,dmmg,ierr)
230: use mex36f90
231: use mex36f90interfaces
233: SNES snes
234: Vec x,f
235: DMMG dmmg
236: PetscErrorCode ierr
238: DMComposite dm
239: Vec fvc1,fvc2
240: PetscInt np,i
241: DA da
242: type(poolfield), pointer :: fHotPool,fColdPool
243: type(channelfield), pointer :: fIHX(:),fCore(:)
244: type(appctx), pointer :: app
245: PetscMPIInt rank
246: type(LocalForm) xlf
248: double precision dhhpl, mult, dhcpl, phpco, pcpio, pcpci, phpii
250: logical debug
252: debug = .false.
254: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
256: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr) ! access the grid information
258: call LocalFormCreate(dm,xlf,ierr) ! Access the local parts of x
259: call LocalFormUpdate(dm,x,xlf,ierr)
261: ! Access the global (non-ghosted) parts of f
262: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
263: call DMCompositeGetAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
264: call DAVecGetArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
265: call DAVecGetArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
267: !BARRY
268: !
269: ! At this point I need to create old time values from "xold" passed in through
270: ! app for
271: !
272: ! xHotPool%vol, xHotPool%vel, xColdPool%vol, xColdPool%vel
273: ! xIHX(i)%press, xCore(i)%press
274: !
275: ! I don't know how to do this.
276: !
277: !BARRY
279: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
280: ! fPool only lives on zeroth process, ghosted xPool lives on all processes
281: if (rank == 0) then
282: !
283: ! compute velocity into hotpool from core
284: !
285: dhhpl = app%dhhpl0 + ( (xlf%HotPool%vol - app%hpvol0) / app%ahp )
286: phpco = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhco) )
287: mult = app%dt / (app%dxc * app%rho)
288: fHotPool%vel = xlf%HotPool%vel - (mult * (xlf%Core(app%nxc-1)%press - phpco) ) + (abs(xlf%HotPool%vel) * xlf%HotPool%vel)
289: !
290: ! compute change in hot pool volume
291: !
292: fHotPool%vol = xlf%HotPool%vol - app%hpvol0 - (app%dt * app%acore * (xlf%HotPool%vel -xlf%ColdPool%vel) )
293: !
294: ! compute velocity into coldpool from IHX
295: !
296: dhcpl = app%dhcpl0 + ( (xlf%ColdPool%vol - app%cpvol0) / app%acp )
297: pcpio = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhio) )
298: mult = app%dt / (app%dxc * app%rho)
299: fColdPool%vel = xlf%ColdPool%vel-(mult*(xlf%IHX(app%nxc-1)%press-pcpio))+(abs(xlf%ColdPool%vel)*xlf%ColdPool%vel)
300: !
301: ! compute the change in cold pool volume
302: !
303: fColdPool%vol = xlf%ColdPool%vol - app%cpvol0 - (app%dt * app%aihx * (xlf%ColdPool%vel - xlf%HotPool%vel) )
304: endif
305: !
306: ! Compute the pressure distribution in IHX and core
307: !
308: pcpci = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhci) )
309: phpii = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhii) )
311: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
313: fIHX(i)%press = xlf%IHX(i)%press - phpii - (app%rho * app%grav * dble(i) * app%dxi)
314: fCore(i)%press = xlf%Core(i)%press - pcpci + (app%rho * app%grav * dble(i) * app%dxc)
316: enddo
318: if (debug) then
319: write(*,*)
320: write(*,*) 'Hot vel,vol ',xlf%HotPool%vel,xlf%HotPool%vol
321: write(*,*) 'delta p = ', xlf%Core(app%nxc-1)%press - phpco,xlf%Core(app%nxc-1)%press,phpco
322: write(*,*)
324: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
325: write(*,*) 'xlf%IHX(',i,')%press = ',xlf%IHX(i)%press
326: enddo
328: write(*,*)
329: write(*,*) 'Cold vel,vol ',xlf%ColdPool%vel,xlf%ColdPool%vol
330: write(*,*) 'delta p = ',xlf%IHX(app%nxc-1)%press - pcpio,xlf%IHX(app%nxc-1)%press, pcpio
331: write(*,*)
334: do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
335: write(*,*) 'xlf%Core(',i,')%press = ',xlf%Core(i)%press
336: enddo
338: endif
340: call DAVecRestoreArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
341: call DAVecRestoreArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
342: call DMCompositeRestoreAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
344: call LocalFormRestore(dm,xlf,ierr)
345: call LocalFormDestroy(dm,xlf,ierr)
346: return
347: end
349: subroutine FormGraph(dmmg,x,view0,view1,ierr)
351: ! --------------------------------------------------------------------------------------
352: !
353: ! FormGraph - Forms Graphics output
354: !
355: ! Input Parameter:
356: ! x - vector
357: ! time - scalar
358: !
359: ! Output Parameters:
360: ! ierr - error code
361: !
362: ! Notes:
363: ! This routine serves as a wrapper for the lower-level routine
364: ! "ApplicationXmgr", where the actual computations are
365: ! done using the standard Fortran style of treating the local
366: ! vector data as a multidimensional array over the local mesh.
367: ! This routine merely accesses the local vector data via
368: ! VecGetArray() and VecRestoreArray().
369: !
370: use mex36f90
371: use mex36f90interfaces
373: Vec x,xvc1,xvc2 !,corep,ihxp
374: DMMG dmmg
375: PetscErrorCode ierr
376: DMComposite dm
377: PetscInt np !,i
378: DA da
379: type(DALocalInfof90) dainfo
380: type(poolfield), pointer :: xHotPool,xColdPool
381: type(channelfield), pointer :: xIHX(:),xCore(:)
382: type(appctx), pointer :: app
384: PetscViewer view0,view1
386: integer iplotnum
387: save iplotnum
388: character*8 grfile
390: data iplotnum / -1 /
391: !BARRY
392: !
393: ! This is my attempt to get the information out of vector "x" to plot
394: ! it. I need to be able to build xCore(i)%press and xIHX(i)%press
395: ! from the vector "x" and I do not know how to do this
396: !
397: !BARRY
399: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
400: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr) ! Access the grid information
402: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
403: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
405: call DMCompositeGetLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
406: call DMCompositeScatter(dm,x,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
407: call DAVecGetArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
408: call DAVecGetArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
410: !
411: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
412: !
413: iplotnum = iplotnum + 1
414: 0
415: !
416: !
417: ! plot corep vector
418: !
419: call VecView(xvc1,view0,ierr);CHKR(ierr)
420: !
421: ! make xmgr plot of corep
422: !
423: write(grfile,4000) iplotnum
424: 4000 format('CoreP',i3.3)
426: open (unit=44,file=grfile,status='unknown')
428: do i = 1,app%nxc
429: write(44,1000) dble(i)*app%dxc, xCore(i)%press
430: enddo
432: close(44)
433: !
434: ! plot ihxp vector
435: !
436: call VecView(xvc2,view1,ierr);CHKR(ierr)
437: !
438: ! make xmgr plot of ihxp
439: !
440: write(grfile,8000) iplotnum
441: 8000 format('IHXPr',i3.3)
443: open (unit=44,file=grfile,status='unknown')
445: do i = 1,app%nxc
446: write(44,1000) dble(i)*app%dxi, xIHX(i)%press
447: enddo
449: close(44)
453: 1000 format(3(e18.12,2x))
455: call DAVecRestoreArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
456: call DAVecRestoreArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
457: call DMCompositeRestoreLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
459: return
460: end
462: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464: subroutine FormInitialGuess(dmmg,v,ierr)
465: use mex36f90
466: use mex36f90interfaces
468: DMMG dmmg
469: Vec v
470: PetscErrorCode ierr
472: DMComposite dm
473: Vec vc1,vc2
474: PetscInt np,i
475: DA da
476: type(DALocalInfof90) dainfo
477: type(poolfield), pointer :: HotPool,ColdPool
478: type(channelfield), pointer :: IHX(:),Core(:)
479: type(appctx), pointer :: app
480: PetscMPIInt rank
482: double precision pcpci, phpii
484: logical debug
486: debug = .false.
488: call DMMGGetUser(dmmg,app,ierr);CHKR(ierr) ! access user context
490: ! Access the grid information
492: call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
493: call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
494: call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)
496: ! Acess the vector unknowns in global (nonghosted) form
498: call DMCompositeGetAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
499: call DAVecGetArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
500: call DAVecGetArrayf90(da,vc2,Core,ierr);CHKR(ierr)
502: call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
503: !
504: ! Set the input values
505: !
507: app%P0 = 1.0d+5
508: app%rho = 1.0d+3
509: app%grav = 9.8d+0
511: app%dhhpl0 = 12.2d+0
512: app%dhcpl0 = 10.16d+0
514: app%lenc = 3.0d+0
515: app%leni = 4.46d+0
517: app%dhci = 0.83d+0
518: app%dhco = app%dhci + app%lenc
520: app%dhii = 7.83d+0
521: app%dhio = app%dhii - app%leni
523: app%dxc = app%lenc / dble(app%nxc)
524: app%dxi = app%leni / dble(app%nxc)
526: app%dt = 1.0d+0
528: app%ahp = 7.0d+0
529: app%acp = 7.0d+0
531: app%acore = 0.8d+0
532: app%aihx = 5.0d-2
534: app%hpvol0 = app%dhhpl0 * app%ahp
535: app%cpvol0 = app%dhcpl0 * app%acp
536: !
537: ! the pools are only unknowns on process 0
538: !
539: if (rank == 0) then
540: HotPool%vel = -1.0d+0
541: HotPool%vol = app%hpvol0
542: ColdPool%vel = 1.0d+0
543: ColdPool%vol = app%cpvol0
544: endif
545: !
546: ! Construct and initial guess at the pressure
547: !
548: pcpci = app%P0 + ( app%rho * app%grav * (app%dhcpl0 - app%dhci) )
549: phpii = app%P0 + ( app%rho * app%grav * (app%dhhpl0 - app%dhii) )
551: if (debug) then
552: write(*,*)
553: write(*,*) 'pcpci = ',pcpci
554: write(*,*) 'phpii = ',phpii
555: write(*,*) 'app%P0 = ',app%P0
556: write(*,*) 'dhcpl0 - app%dhci ',app%dhcpl0 - app%dhci,app%dhcpl0, app%dhci
557: write(*,*) 'dhhpl0 - app%dhii ',app%dhhpl0 - app%dhii,app%dhhpl0, app%dhii
558: write(*,*)
559: endif
561: do i=dainfo%xs,dainfo%xs+dainfo%xm-1
563: IHX(i)%press = phpii + (app%rho * app%grav * dble(i) * app%dxi)
565: Core(i)%press = pcpci - (app%rho * app%grav * dble(i) * app%dxc)
567: enddo
569: call DAVecRestoreArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
570: call DAVecRestoreArrayf90(da,vc2,Core,ierr);CHKR(ierr)
571: call DMCompositeRestoreAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
573: CHKMEMQ
574: return
575: end