Actual source code: ex36f90.F90

  1: #define PETSC_AVOID_DECLARATIONS
 2:  #include include/finclude/petscall.h

  4: !  Error handler that aborts when error is detected
  5: !
  6:       subroutine HE1(ierr,line)
  7:       use mex36f90
  8:       use mex36f90interfaces

 10:       call PetscError(ierr,line,0,'',ierr)
 11:       call MPI_Abort(PETSC_COMM_WORLD,ierr,ierr)
 12:       return
 13:       end
 14: #define CHKQ(n) if(n .ne. 0)call HE1(n,__LINE__)

 16: !  Error handler forces traceback of where error occurred
 17: !
 18:       subroutine HE2(ierr,line)
 19:       use mex36f90
 20:       use mex36f90interfaces

 22:       call PetscError(ierr,line,0,'',ierr)
 23:       return
 24:       end
 25: #define CHKR(n) if(n .ne. 0)then;call HE2(n,__LINE__);return;endif

 27: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 28: !   Utilities to map between global (linear algebra) to local (ghosted, physics) representations
 29: !
 30: !    Allocates the local work arrays and vectors
 31:       subroutine LocalFormCreate(dm,lf,ierr)
 32:       use mex36f90
 33:       use mex36f90interfaces
 34:       implicit none
 35:       DMComposite  dm
 36:       type(LocalForm) lf
 37:       PetscErrorCode ierr

 39:       call DMCompositeGetEntries(dm,lf%da,lf%np,lf%da,lf%np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
 40:       call DAGetLocalInfof90(lf%da,lf%dainfo,ierr);CHKR(ierr)

 42:       call DMCompositeGetLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 43:       CHKR(ierr)
 44:       return
 45:       end

 47: !     Updates the (ghosted) IHX and Core arrays from the global vector x
 48:       subroutine LocalFormUpdate(dm,x,lf,ierr)
 49:       use mex36f90
 50:       use mex36f90interfaces
 51:       implicit none
 52:       DMComposite  dm
 53:       type(LocalForm) lf
 54:       PetscErrorCode ierr
 55:       Vec x

 57:       call DMCompositeScatter(dm,x,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 58:       CHKR(ierr)
 59:       call DAVecGetArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
 60:       call DAVecGetArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
 61:       return
 62:       end

 64: !     You should call this when you no longer need access to %IHX and %Core
 65:       subroutine LocalFormRestore(dm,lf,ierr)
 66:       use mex36f90
 67:       use mex36f90interfaces
 68:       implicit none
 69:       DMComposite  dm
 70:       type(LocalForm) lf
 71:       PetscErrorCode ierr

 73:       call DAVecRestoreArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
 74:       call DAVecRestoreArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
 75:       return
 76:       end

 78: !     Destroys the data structure
 79:       subroutine LocalFormDestroy(dm,lf,ierr)
 80:       use mex36f90
 81:       use mex36f90interfaces
 82:       implicit none
 83:       DMComposite  dm
 84:       type(LocalForm) lf
 85:       PetscErrorCode ierr

 87:       call DMCompositeRestoreLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 88:       CHKR(ierr)
 89:       return
 90:       end

 92: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 93: !      Implements a nonlinear solver for a simple domain
 94: !     consisting of 2 zero dimensional problems linked by
 95: !     2 one dimensional problems.
 96: !
 97: !                           Channel1
 98: !                       -------------------------
 99: !               Pool 1                              Pool 2
100: !                       -------------------------
101: !                           Channel2
102: !VAM
103: !VAM
104: !VAM
105: !VAM                         Hot Pool 1
106: !VAM                 |                       |
107: !VAM                 |                       |
108: !VAM                 |                       |
109: !VAM                 |                       |
110: !VAM  Core Channel 2 |                       | IHX Channel 1
111: !VAM                 |                       |
112: !VAM                 |                       |
113: !VAM                 |                       |
114: !VAM                 |                       |
115: !VAM                         Cold Pool 2
116: !VAM
117: !
118: !     For Newton's method with block Jacobi (one block for
119: !     each channel and one block for each pool) run with
120: !
121: !       -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
122: !       -help lists all options

124:       program ex36f90
125:       use mex36f90
126: !     use mex36f90interfaces
127:  #include include/finclude/petsc.h
128:  #include include/finclude/petscviewer.h
129:  #include include/finclude/petscvec.h


132:       DMMGArray        dmmg
133:       DMMG             dmmglevel
134:       PetscErrorCode   ierr
135:       DA               da
136:       DMComposite      dm
137:       type(appctx)     app
138:       external         FormInitialGuess,FormFunction,FormGraph
139:       Vec              x

141:       double precision time
142:       integer i
143: !BARRY
144:        PetscViewer        view0,view1
145: !BARRY

147:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

149: !      Create the composite DM object that manages the grid

151:       call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
152:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
153:       CHKR(ierr)
154:       call DMCompositeAddDA(dm,da,ierr);CHKR(ierr)
155:       call DADestroy(da,ierr);CHKR(ierr)
156:       call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
157:       call DMCompositeAddDA(dm,da,ierr);CHKR(ierr)
158:       call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)

160: !       Create solver object and associate it with the unknowns (on the grid)

162:       call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
163:       call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
164:       call DMCompositeDestroy(dm,ierr);CHKR(ierr)
165:       call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr)  ! currently only one level solver
166:       call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
167:       call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(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