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