Actual source code: ex6f90.F

petsc-3.4.5 2014-06-29
  1: !-----------------------------------------------------------------------
  2: !
  3: !     manages tokamak edge region.
  4: !
  5: !     This is a translation of ex6.c into F90 by Alan Glasser, LANL
  6: !-----------------------------------------------------------------------
  7: !-----------------------------------------------------------------------
  8: !     code organization.
  9: !-----------------------------------------------------------------------
 10: !     0. barry_mod.
 11: !     1. barry_get_global_corners.
 12: !     2. barry_get_global_vector.
 13: !     3. barry_get_local_vector.
 14: !     4. barry_global_to_local.
 15: !     5. barry_destroy_fa.
 16: !     6. barry_create_fa.
 17: !     7. barry_draw_patch.
 18: !     8. barry_draw_fa.
 19: !     9. barry_map_region3.
 20: !     10. barry_map_region2.
 21: !     11. barry_map_region1.
 22: !     12. barry_main.
 23: !-----------------------------------------------------------------------
 24: !     subprogram 0. barry_mod.
 25: !     module declarations.
 26: !-----------------------------------------------------------------------
 27: !-----------------------------------------------------------------------
 28: !     declarations.
 29: !-----------------------------------------------------------------------
 30:       MODULE barry_mod
 31:       IMPLICIT NONE

 33: #include <finclude/petscsys.h>
 34: #include <finclude/petscvec.h>
 35: #include <finclude/petscdmda.h>
 36: #include <finclude/petscmat.h>
 37: #include <finclude/petscis.h>
 38: #include <finclude/petscao.h>
 39: #include <finclude/petscviewer.h>
 40: #include <finclude/petscdraw.h>

 42:       TYPE :: fa_type
 43:       MPI_Comm, DIMENSION(0:2) :: comm
 44:       INTEGER, DIMENSION(0:2) :: xl,yl,ml,nl,xg,yg,mg,ng,offl,offg
 45:       Vec :: g,l
 46:       VecScatter :: vscat
 47:       INTEGER :: p1,p2,r1,r2,r1g,r2g,sw
 48:       END TYPE fa_type

 50:       TYPE :: patch_type
 51:       INTEGER :: mx,my
 52:       PetscScalar, DIMENSION(:,:,:), POINTER :: xy
 53:       END TYPE patch_type

 55:       LOGICAL, PARAMETER :: diagnose=.FALSE.
 56:       INTEGER :: ierr
 57:       REAL(8), PARAMETER :: pi=3.1415926535897932385_8,twopi=2*pi

 59:       CONTAINS
 60: !-----------------------------------------------------------------------
 61: !     subprogram 1. barry_get_global_corners.
 62: !     returns global corner data.
 63: !-----------------------------------------------------------------------
 64: !-----------------------------------------------------------------------
 65: !     declarations.
 66: !-----------------------------------------------------------------------
 67:       SUBROUTINE barry_get_global_corners(fa,j,x,y,m,n)

 69:       TYPE(fa_type), INTENT(IN) :: fa
 70:       INTEGER, INTENT(IN) :: j
 71:       INTEGER, INTENT(OUT) :: x,y,m,n
 72: !-----------------------------------------------------------------------
 73: !     computations.
 74: !-----------------------------------------------------------------------
 75:       IF (fa%comm(j) /= 0) THEN
 76:          x=fa%xg(j)
 77:          y=fa%yg(j)
 78:          m=fa%mg(j)
 79:          n=fa%ng(j)
 80:       ELSE
 81:          x=0
 82:          y=0
 83:          m=0
 84:          n=0
 85:       ENDIF
 86: !-----------------------------------------------------------------------
 87: !     terminate.
 88: !-----------------------------------------------------------------------
 89:       RETURN
 90:       END SUBROUTINE barry_get_global_corners
 91: !-----------------------------------------------------------------------
 92: !     subprogram 2. barry_get_global_vector.
 93: !     returns local vector.
 94: !-----------------------------------------------------------------------
 95: !-----------------------------------------------------------------------
 96: !     declarations.
 97: !-----------------------------------------------------------------------
 98:       SUBROUTINE barry_get_global_vector(fa,v)

100:       TYPE(fa_type), INTENT(IN) :: fa
101:       Vec, INTENT(OUT) :: v
102: !-----------------------------------------------------------------------
103: !     computations.
104: !-----------------------------------------------------------------------
105:       CALL VecDuplicate(fa%g,v,ierr)
106: !-----------------------------------------------------------------------
107: !     terminate.
108: !-----------------------------------------------------------------------
109:       RETURN
110:       END SUBROUTINE barry_get_global_vector
111: !-----------------------------------------------------------------------
112: !     subprogram 3. barry_get_local_vector.
113: !     returns local vector.
114: !-----------------------------------------------------------------------
115: !-----------------------------------------------------------------------
116: !     declarations.
117: !-----------------------------------------------------------------------
118:       SUBROUTINE barry_get_local_vector(fa,v)

120:       TYPE(fa_type), INTENT(IN) :: fa
121:       Vec, INTENT(OUT) :: v
122: !-----------------------------------------------------------------------
123: !     computations.
124: !-----------------------------------------------------------------------
125:       CALL VecDuplicate(fa%l,v,ierr)
126: !-----------------------------------------------------------------------
127: !     terminate.
128: !-----------------------------------------------------------------------
129:       RETURN
130:       END SUBROUTINE barry_get_local_vector
131: !-----------------------------------------------------------------------
132: !     subprogram 4. barry_global_to_local.
133: !     performs VecScatter.
134: !-----------------------------------------------------------------------
135: !-----------------------------------------------------------------------
136: !     declarations.
137: !-----------------------------------------------------------------------
138:       SUBROUTINE barry_global_to_local(fa,g,l)

140:       TYPE(fa_type), INTENT(IN) :: fa
141:       Vec, INTENT(IN) :: g
142:       Vec, INTENT(OUT) :: l
143: !-----------------------------------------------------------------------
144: !     computations.
145: !-----------------------------------------------------------------------
146:       CALL VecScatterBegin(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD,     &
147:      &     ierr)
148:       CALL VecScatterEnd(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD,       &
149:      &     ierr)
150: !-----------------------------------------------------------------------
151: !     terminate.
152: !-----------------------------------------------------------------------
153:       RETURN
154:       END SUBROUTINE barry_global_to_local
155: !-----------------------------------------------------------------------
156: !     subprogram 5. barry_destroy_fa.
157: !     destroys fa.
158: !-----------------------------------------------------------------------
159: !-----------------------------------------------------------------------
160: !     declarations.
161: !-----------------------------------------------------------------------
162:       SUBROUTINE barry_destroy_fa(fa)

164:       TYPE(fa_type), INTENT(OUT) :: fa
165: !-----------------------------------------------------------------------
166: !     computations.
167: !-----------------------------------------------------------------------
168:       CALL VecDestroy(fa%g,ierr)
169:       CALL VecDestroy(fa%l,ierr)
170:       CALL VecScatterDestroy(fa%vscat,ierr)
171: !-----------------------------------------------------------------------
172: !     terminate.
173: !-----------------------------------------------------------------------
174:       RETURN
175:       END SUBROUTINE barry_destroy_fa
176: !-----------------------------------------------------------------------
177: !     subprogram 6. barry_create_fa.
178: !     creates fa.
179: !-----------------------------------------------------------------------
180: !-----------------------------------------------------------------------
181: !     declarations.
182: !-----------------------------------------------------------------------
183:       SUBROUTINE barry_create_fa(fa)

185:       TYPE(fa_type), INTENT(OUT) :: fa

187:       INTEGER :: rank,tonglobal,globalrstart,x,nx,y,ny,i,j,k,ig,              &
188:      &     fromnglobal,nscat,nlocal,cntl1,cntl2,cntl3,il,it
189:       INTEGER, DIMENSION(0:2) :: offt
190:       INTEGER, DIMENSION(:), POINTER :: tonatural,fromnatural,                &
191:      &     to,from,indices
192:       PetscScalar, DIMENSION(1) :: globalarray
193:       PetscScalar, DIMENSION(1) :: localarray
194:       PetscScalar, DIMENSION(1) :: toarray

196:       DM :: da1=0,da2=0,da3=0
197:       Vec :: vl1=0,vl2=0,vl3=0
198:       Vec :: vg1=0,vg2=0,vg3=0
199:       AO :: toao,globalao
200:       IS :: tois,globalis,is
201:       Vec :: tovec,globalvec,localvec
202:       VecScatter :: vscat
203: !-----------------------------------------------------------------------
204: !     set integer members of fa.
205: !-----------------------------------------------------------------------
206:       fa%p1=24
207:       fa%p2=15
208:       fa%r1=6
209:       fa%r2=6
210:       fa%sw=1
211:       fa%r1g=fa%r1+fa%sw
212:       fa%r2g=fa%r2+fa%sw
213: !-----------------------------------------------------------------------
214: !     set up communicators.
215: !-----------------------------------------------------------------------
216:       CALL MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
217:       fa%comm=PETSC_COMM_WORLD
218: !-----------------------------------------------------------------------
219: !     set up distributed arrays.
220: !-----------------------------------------------------------------------
221:       IF (fa%comm(1) /= 0) THEN
222:          CALL DMDACreate2d(fa%comm(1),DMDA_BOUNDARY_PERIODIC,                   &
223:      &        DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,                          &
224:      &        fa%p2,fa%r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw,                &
225:      &        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da2,ierr)
226:          CALL DMGetLocalVector(da2,vl2,ierr)
227:          CALL DMGetGlobalVector(da2,vg2,ierr)
228:       ENDIF
229:       IF (fa%comm(2) /= 0) THEN
230:          CALL DMDACreate2d(fa%comm(2),DMDA_BOUNDARY_NONE,                       &
231:      &        DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,                        &
232:      &        fa%p1-fa%p2,fa%r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw,          &
233:      &        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da3,ierr)
234:          CALL DMGetLocalVector(da3,vl3,ierr)
235:          CALL DMGetGlobalVector(da3,vg3,ierr)
236:       ENDIF
237:       IF (fa%comm(0) /= 0) THEN
238:          CALL DMDACreate2d(fa%comm(0),DMDA_BOUNDARY_NONE,                       &
239:      &        DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,                        &
240:      &        fa%p1,fa%r1g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw,                &
241:      &        PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da1,ierr)
242:          CALL DMGetLocalVector(da1,vl1,ierr)
243:          CALL DMGetGlobalVector(da1,vg1,ierr)
244:       ENDIF
245: !-----------------------------------------------------------------------
246: !     count the number of unknowns owned on each processor and determine
247: !     the starting point of each processors ownership
248: !     for global vector with redundancy.
249: !-----------------------------------------------------------------------
250:       tonglobal = 0
251:       IF (fa%comm(1) /= 0) THEN
252:          CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
253:          tonglobal=tonglobal+nx*ny
254:       ENDIF
255:       IF (fa%comm(2) /= 0) THEN
256:          CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
257:          tonglobal=tonglobal+nx*ny
258:       ENDIF
259:       IF (fa%comm(0) /= 0) THEN
260:          CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
261:          tonglobal=tonglobal+nx*ny
262:       ENDIF
263:       WRITE(*,'("[",i1,"]",a,i3)')                                               &
264:      &     rank," Number of unknowns owned ",tonglobal
265: !-----------------------------------------------------------------------
266: !     get tonatural number for each node
267: !-----------------------------------------------------------------------
268:       ALLOCATE(tonatural(0:tonglobal))
269:       tonglobal=0
270:       IF (fa%comm(1) /= 0) THEN
271:          CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
272:          DO j=0,ny-1
273:             DO i=0,nx-1
274:                tonatural(tonglobal)=(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
275:                tonglobal=tonglobal+1
276:             ENDDO
277:          ENDDO
278:       ENDIF
279:       IF (fa%comm(2) /= 0) THEN
280:          CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
281:          DO j=0,ny-1
282:             DO i=0,nx-1
283:                IF (x+i < (fa%p1-fa%p2)/2) THEN
284:                   tonatural(tonglobal)=x+i+fa%p1*(y+j)
285:                ELSE
286:                   tonatural(tonglobal)=fa%p2+x+i+fa%p1*(y+j)
287:                ENDIF
288:                tonglobal=tonglobal+1
289:             ENDDO
290:          ENDDO
291:       ENDIF
292:       IF (fa%comm(0) /= 0) THEN
293:          CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
294:          DO j=0,ny-1
295:             DO i=0,nx-1
296:                tonatural(tonglobal)=fa%p1*fa%r2g+x+i+fa%p1*(y+j)
297:                tonglobal=tonglobal+1
298:             ENDDO
299:          ENDDO
300:       ENDIF
301: !-----------------------------------------------------------------------
302: !     diagnose tonatural.
303: !-----------------------------------------------------------------------
304:       IF (diagnose)THEN
305:          WRITE(*,'(a,i3,a)')"tonglobal = ",tonglobal,", tonatural = "
306:          WRITE(*,'(2i4)')(i,tonatural(i),i=0,tonglobal-1)
307:       ENDIF
308: !-----------------------------------------------------------------------
309: !     create application ordering and deallocate tonatural.
310: !-----------------------------------------------------------------------
311:       CALL AOCreateBasic(PETSC_COMM_WORLD,tonglobal,tonatural,                   &
312:      &     PETSC_NULL_INTEGER,toao,ierr)
313:       DEALLOCATE(tonatural)
314: !-----------------------------------------------------------------------
315: !     count the number of unknowns owned on each processor and determine
316: !     the starting point of each processors ownership for global vector
317: !     without redundancy.
318: !-----------------------------------------------------------------------
319:       fromnglobal=0
320:       fa%offg(1)=0
321:       offt(1)=0
322:       IF (fa%comm(1) /= 0)THEN
323:          CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
324:          offt(2)=nx*ny
325:          IF (y+ny == fa%r2g)ny=ny-1
326:          fromnglobal=fromnglobal+nx*ny
327:          fa%offg(2)=fromnglobal
328:       ELSE
329:          offt(2)=0
330:          fa%offg(2)=0
331:       ENDIF
332:       IF (fa%comm(2) /= 0)THEN
333:          CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
334:          offt(0)=offt(2)+nx*ny
335:          IF (y+ny == fa%r2g)ny=ny-1
336:          fromnglobal=fromnglobal+nx*ny
337:          fa%offg(0)=fromnglobal
338:       ELSE
339:          offt(0)=offt(2)
340:          fa%offg(0)=fromnglobal
341:       ENDIF
342:       IF (fa%comm(0) /= 0)THEN
343:          CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
344:          IF (y == 0)ny=ny-1
345:          fromnglobal=fromnglobal+nx*ny
346:       ENDIF
347:       CALL MPI_Scan(fromnglobal,globalrstart,1,MPI_INTEGER,                      &
348:      &     MPI_SUM,PETSC_COMM_WORLD,ierr)
349:       globalrstart=globalrstart-fromnglobal
350:       WRITE(*,'("[",i1,"]",a,i3)')                                               &
351:      &     rank," Number of unknowns owned ",fromnglobal
352: !-----------------------------------------------------------------------
353: !     get fromnatural number for each node.
354: !-----------------------------------------------------------------------
355:       ALLOCATE(fromnatural(0:fromnglobal))
356:       fromnglobal=0
357:       IF (fa%comm(1) /= 0)THEN
358:          CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
359:          IF (y+ny == fa%r2g)ny=ny-1
360:          fa%xg(1)=x
361:          fa%yg(1)=y
362:          fa%mg(1)=nx
363:          fa%ng(1)=ny
364:          CALL DMDAGetGhostCorners(da2,fa%xl(1),fa%yl(1),0,fa%ml(1),                &
365:      &        fa%nl(1),0,ierr)
366:          DO j=0,ny-1
367:             DO i=0,nx-1
368:                fromnatural(fromnglobal)                                          &
369:      &              =(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
370:                fromnglobal=fromnglobal+1
371:             ENDDO
372:          ENDDO
373:       ENDIF
374:       IF (fa%comm(2) /= 0)THEN
375:          CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
376:          IF (y+ny == fa%r2g)ny=ny-1
377:          fa%xg(2)=x
378:          fa%yg(2)=y
379:          fa%mg(2)=nx
380:          fa%ng(2)=ny
381:          CALL DMDAGetGhostCorners(da3,fa%xl(2),fa%yl(2),0,fa%ml(2),                &
382:      &        fa%nl(2),0,ierr)
383:          DO j=0,ny-1
384:             DO i=0,nx-1
385:                IF (x+i < (fa%p1-fa%p2)/2)THEN
386:                   fromnatural(fromnglobal)=x+i+fa%p1*(y+j)
387:                ELSE
388:                   fromnatural(fromnglobal)=fa%p2+x+i+fa%p1*(y+j)
389:                ENDIF
390:                fromnglobal=fromnglobal+1
391:             ENDDO
392:          ENDDO
393:       ENDIF
394:       IF (fa%comm(0) /= 0)THEN
395:          CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
396:          IF (y == 0)THEN
397:             ny=ny-1
398:          ELSE
399:             y=y-1
400:          ENDIF
401:          fa%xg(0)=x
402:          fa%yg(0)=y
403:          fa%mg(0)=nx
404:          fa%ng(0)=ny
405:          CALL DMDAGetGhostCorners(da1,fa%xl(0),fa%yl(0),0,fa%ml(0),               &
406:      &        fa%nl(0),0,ierr)
407:          DO j=0,ny-1
408:             DO i=0,nx-1
409:                fromnatural(fromnglobal)=fa%p1*fa%r2+x+i+fa%p1*(y+j)
410:                fromnglobal=fromnglobal+1
411:             ENDDO
412:          ENDDO
413:       ENDIF
414: !-----------------------------------------------------------------------
415: !     diagnose fromnatural.
416: !-----------------------------------------------------------------------
417:       IF (diagnose)THEN
418:          WRITE(*,'(a,i3,a)')"fromnglobal = ",fromnglobal                        &
419:      &        ,", fromnatural = "
420:          WRITE(*,'(2i4)')(i,fromnatural(i),i=0,fromnglobal-1)
421:       ENDIF
422: !-----------------------------------------------------------------------
423: !     create application ordering and deallocate fromnatural.
424: !-----------------------------------------------------------------------
425:       CALL AOCreateBasic(PETSC_COMM_WORLD,fromnglobal,fromnatural,              &
426:      &     PETSC_NULL_INTEGER,globalao,ierr)
427:       DEALLOCATE(fromnatural)
428: !-----------------------------------------------------------------------
429: !     set up scatter that updates 1 from 2 and 3 and 3 and 2 from 1
430: !-----------------------------------------------------------------------
431:       ALLOCATE(to(0:tonglobal),from(0:tonglobal))
432:       nscat=0
433:       IF (fa%comm(1) /= 0)THEN
434:          CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
435:          DO j=0,ny-1
436:             DO i=0,nx-1
437:                to(nscat)=(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
438:                from(nscat)=to(nscat)
439:                nscat=nscat+1
440:             ENDDO
441:          ENDDO
442:       ENDIF
443:       IF (fa%comm(2) /= 0)THEN
444:          CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
445:          DO j=0,ny-1
446:             DO i=0,nx-1
447:                IF (x+i < (fa%p1-fa%p2)/2)THEN
448:                   to(nscat)=x+i+fa%p1*(y+j)
449:                ELSE
450:                   to(nscat)=fa%p2+x+i+fa%p1*(y+j)
451:                ENDIF
452:                from(nscat)=to(nscat)
453:                nscat=nscat+1
454:             ENDDO
455:          ENDDO
456:       ENDIF
457:       IF (fa%comm(0) /= 0)THEN
458:          CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
459:          DO j=0,ny-1
460:             DO i=0,nx-1
461:                to(nscat)=fa%p1*fa%r2g+x+i+fa%p1*(y+j)
462:                from(nscat)=fa%p1*(fa%r2-1)+x+i+fa%p1*(y+j)
463:                nscat=nscat+1
464:             ENDDO
465:          ENDDO
466:       ENDIF
467: !-----------------------------------------------------------------------
468: !     diagnose to and from.
469: !-----------------------------------------------------------------------
470:       IF (diagnose)THEN
471:          WRITE(*,'(a,i3,a)')"nscat = ",nscat,", to, from = "
472:          WRITE(*,'(3i4)')(i,to(i),from(i),i=0,nscat-1)
473:       ENDIF
474: !-----------------------------------------------------------------------
475: !     create vecscatter.
476: !-----------------------------------------------------------------------
477:       CALL AOApplicationToPetsc(toao,nscat,to,ierr)
478:       CALL AOApplicationToPetsc(globalao,nscat,from,ierr)
479:       CALL ISCreateGeneral(PETSC_COMM_WORLD,nscat,to,PETSC_COPY_VALUES,               &
480:      &                     tois,ierr)
481:       CALL ISCreateGeneral(PETSC_COMM_WORLD,nscat,from,PETSC_COPY_VALUES              &
482:      &                     ,globalis,ierr)
483:       CALL VecCreateMPI(PETSC_COMM_WORLD,tonglobal,PETSC_DETERMINE,                   &
484:      &     tovec,ierr)
485:       CALL VecCreateMPI(PETSC_COMM_WORLD,fromnglobal,PETSC_DETERMINE,                 &
486:      &     globalvec,ierr)
487:       CALL VecScatterCreate(globalvec,globalis,tovec,tois,vscat,ierr)
488: !-----------------------------------------------------------------------
489: !     clean up.
490: !-----------------------------------------------------------------------
491:       CALL ISDestroy(tois,ierr)
492:       CALL ISDestroy(globalis,ierr)
493:       CALL AODestroy(globalao,ierr)
494:       CALL AODestroy(toao,ierr)
495:       DEALLOCATE(to,from)
496: !-----------------------------------------------------------------------
497: !     fill up global vector without redundant values with PETSc global
498: !     numbering
499: !-----------------------------------------------------------------------
500:       CALL VecGetArray(globalvec,globalarray,ig,ierr)
501:       k=ig
502:       IF (diagnose)WRITE(*,'(a,i3,a)')                                                &
503:      &     "fromnglobal = ",fromnglobal,", globalarray = "
504:       DO i=0,fromnglobal-1
505:          k=k+1
506:          globalarray(k)=globalrstart+i
507:          IF (diagnose)WRITE(*,'(i4,f11.3)')i,globalarray(k)
508:       ENDDO
509:       CALL VecRestoreArray(globalvec,globalarray,ig,ierr)
510: !-----------------------------------------------------------------------
511: !     scatter PETSc global indices to redundant valued array.
512: !-----------------------------------------------------------------------
513:       CALL VecScatterBegin(vscat,globalvec,tovec,INSERT_VALUES,                      &
514:      &     SCATTER_FORWARD,ierr)
515:       CALL VecScatterEnd(vscat,globalvec,tovec,INSERT_VALUES,                        &
516:      &     SCATTER_FORWARD,ierr)
517: !-----------------------------------------------------------------------
518: !     create local vector as concatenation of local vectors
519: !-----------------------------------------------------------------------
520:       nlocal=0
521:       cntl1=0
522:       cntl2=0
523:       cntl3=0
524:       IF (fa%comm(1) /= 0)THEN
525:          CALL VecGetSize(vl2,cntl2,ierr)
526:          nlocal=nlocal+cntl2
527:       ENDIF
528:       IF (fa%comm(2) /= 0)THEN
529:          CALL VecGetSize(vl3,cntl3,ierr)
530:          nlocal=nlocal+cntl3
531:       ENDIF
532:       IF (fa%comm(0) /= 0)THEN
533:          CALL VecGetSize(vl1,cntl1,ierr)
534:          nlocal=nlocal+cntl1
535:       ENDIF
536:       fa%offl(0)=cntl2+cntl3
537:       fa%offl(1)=0
538:       fa%offl(2)=cntl2
539:       CALL VecCreateSeq(PETSC_COMM_SELF,nlocal,localvec,ierr)
540: !-----------------------------------------------------------------------
541: !     cheat so that vl1, vl2, vl3 shared array memory with localvec.
542: !-----------------------------------------------------------------------
543:       CALL VecGetArray(localvec,localarray,il,ierr)
544:       CALL VecGetArray(tovec,toarray,it,ierr)
545:       IF (fa%comm(1) /= 0)THEN
546:          CALL VecPlaceArray(vl2,localarray(il+1+fa%offl(1)),ierr)
547:          CALL VecPlaceArray(vg2,toarray(it+1+offt(1)),ierr)
548:          CALL DMGlobalToLocalBegin(da2,vg2,INSERT_VALUES,vl2,ierr)
549:          CALL DMGlobalToLocalEnd(da2,vg2,INSERT_VALUES,vl2,ierr)
550:          CALL DMRestoreGlobalVector(da2,vg2,ierr)
551:       ENDIF
552:       IF (fa%comm(2) /= 0)THEN
553:          CALL VecPlaceArray(vl3,localarray(il+1+fa%offl(2)),ierr)
554:          CALL VecPlaceArray(vg3,toarray(it+1+offt(2)),ierr)
555:          CALL DMGlobalToLocalBegin(da3,vg3,INSERT_VALUES,vl3,ierr)
556:          CALL DMGlobalToLocalEnd(da3,vg3,INSERT_VALUES,vl3,ierr)
557:          CALL DMRestoreGlobalVector(da3,vg3,ierr)
558:       ENDIF
559:       IF (fa%comm(0) /= 0)THEN
560:          CALL VecPlaceArray(vl1,localarray(il+1+fa%offl(0)),ierr)
561:          CALL VecPlaceArray(vg1,toarray(it+1+offt(0)),ierr)
562:          CALL DMGlobalToLocalBegin(da1,vg1,INSERT_VALUES,vl1,ierr)
563:          CALL DMGlobalToLocalEnd(da1,vg1,INSERT_VALUES,vl1,ierr)
564:          CALL DMRestoreGlobalVector(da1,vg1,ierr)
565:       ENDIF
566:       CALL VecRestoreArray(localvec,localarray,il,ierr)
567:       CALL VecRestoreArray(tovec,toarray,it,ierr)
568: !-----------------------------------------------------------------------
569: !     clean up.
570: !-----------------------------------------------------------------------
571:       CALL VecScatterDestroy(vscat,ierr)
572:       CALL VecDestroy(tovec,ierr)
573: !-----------------------------------------------------------------------
574: !     compute index set.
575: !-----------------------------------------------------------------------
576:       ALLOCATE(indices(0:nlocal-1))
577:       CALL VecGetArray(localvec,localarray,il,ierr)
578:       k=il
579:       IF (diagnose)WRITE(*,'(a,i3,a3,i4,a)')                                     &
580:      &     "nlocal = ", nlocal,", offl = ",fa%offl,", indices = "
581:       DO i=0,nlocal-1
582:          k=k+1
583:          indices(i)= int(localarray(k))
584:          IF (diagnose)WRITE(*,'(2i4)')i,indices(i)
585:       ENDDO
586:       CALL VecRestoreArray(localvec,localarray,il,ierr)
587:       CALL ISCreateBlock(PETSC_COMM_WORLD,2,nlocal,indices,                     &
588:      &                   PETSC_COPY_VALUES,is,ierr)
589:       DEALLOCATE(indices)
590: !-----------------------------------------------------------------------
591: !     create local and global vectors.
592: !-----------------------------------------------------------------------
593:       CALL VecCreateSeq(PETSC_COMM_SELF,2*nlocal,fa%l,ierr)
594:       CALL VecCreateMPI(PETSC_COMM_WORLD,2*fromnglobal,PETSC_DETERMINE,         &
595:      &     fa%g,ierr)
596: !-----------------------------------------------------------------------
597: !     create final scatter that goes directly from globalvec to localvec
598: !     this is the one to be used in the application code.
599: !-----------------------------------------------------------------------
600:       CALL VecScatterCreate(fa%g,is,fa%l,PETSC_NULL_OBJECT,                     &
601:      &     fa%vscat,ierr)
602: !-----------------------------------------------------------------------
603: !     clean up.
604: !-----------------------------------------------------------------------
605:       CALL ISDestroy(is,ierr)
606:       CALL VecDestroy(globalvec,ierr)
607:       CALL VecDestroy(localvec,ierr)
608:       IF (fa%comm(0) /= 0)THEN
609:          CALL DMRestoreLocalVector(da1,vl1,ierr)
610:          CALL DMDestroy(da1,ierr)
611:       ENDIF
612:       IF (fa%comm(1) /= 0)THEN
613:          CALL DMRestoreLocalVector(da2,vl2,ierr)
614:          CALL DMDestroy(da2,ierr)
615:       ENDIF
616:       IF (fa%comm(2) /= 0)THEN
617:          CALL DMRestoreLocalVector(da3,vl3,ierr)
618:          CALL DMDestroy(da3,ierr)
619:       ENDIF
620: !-----------------------------------------------------------------------
621: !     terminate.
622: !-----------------------------------------------------------------------
623:       RETURN
624:       END SUBROUTINE barry_create_fa
625: !-----------------------------------------------------------------------
626: !     subprogram 7. barry_draw_patch.
627: !     crude graphics to test that the ghost points are properly updated.
628: !-----------------------------------------------------------------------
629: !-----------------------------------------------------------------------
630: !     declarations.
631: !-----------------------------------------------------------------------
632:       SUBROUTINE barry_draw_patch(draw,patch,ierr)

634:       PetscDraw, INTENT(IN) :: draw
635:       TYPE(patch_type), DIMENSION(0:2), INTENT(IN) :: patch
636:       INTEGER, INTENT(OUT) :: ierr

638:       INTEGER :: ix,iy,j
639:       PetscReal, DIMENSION(4) :: x,y
640: !-----------------------------------------------------------------------
641: !     draw it.
642: !-----------------------------------------------------------------------
643:       DO j=0,2
644:          DO iy=1,patch(j)%my
645:             DO ix=1,patch(j)%mx
646:                x(1)=patch(j)%xy(1,ix-1,iy-1)
647:                y(1)=patch(j)%xy(2,ix-1,iy-1)
648:                x(2)=patch(j)%xy(1,ix,iy-1)
649:                y(2)=patch(j)%xy(2,ix,iy-1)
650:                x(3)=patch(j)%xy(1,ix,iy)
651:                y(3)=patch(j)%xy(2,ix,iy)
652:                x(4)=patch(j)%xy(1,ix-1,iy)
653:                y(4)=patch(j)%xy(2,ix-1,iy)
654:                CALL PetscDrawLine(draw,x(1),y(1),x(2),y(2),                        &
655:      &              PETSC_DRAW_BLACK,ierr)
656:                CALL PetscDrawLine(draw,x(2),y(2),x(3),y(3),                        &
657:      &              PETSC_DRAW_BLACK,ierr)
658:                CALL PetscDrawLine(draw,x(3),y(3),x(4),y(4),                        &
659:      &              PETSC_DRAW_BLACK,ierr)
660:                CALL PetscDrawLine(draw,x(4),y(4),x(1),y(1),                        &
661:      &              PETSC_DRAW_BLACK,ierr)
662:             ENDDO
663:          ENDDO
664:       ENDDO
665: !-----------------------------------------------------------------------
666: !     terminate.
667: !-----------------------------------------------------------------------
668:       ierr=0
669:       RETURN
670:       END SUBROUTINE barry_draw_patch
671: !-----------------------------------------------------------------------
672: !     subprogram 8. barry_draw_fa.
673: !     deallocates local array.
674: !-----------------------------------------------------------------------
675: !-----------------------------------------------------------------------
676: !     declarations.
677: !-----------------------------------------------------------------------
678:       SUBROUTINE barry_draw_fa(fa,v)

680:       TYPE(fa_type), INTENT(IN) :: fa
681:       Vec, INTENT(IN) :: v

683:       INTEGER :: iv,vn,ln,j,offset
684:       REAL(8), DIMENSION(1) :: va
685:       TYPE(patch_type), DIMENSION(0:2) :: patch
686:       PetscDraw :: draw
687:       PetscReal :: xmin,xmax,ymin,ymax
688:       PetscReal :: ymint,xmint,ymaxt,xmaxt
689: !-----------------------------------------------------------------------
690: !     set extrema.
691: !-----------------------------------------------------------------------
692:       xmin=HUGE(xmin)
693:       xmax=-HUGE(xmax)
694:       ymin=HUGE(ymin)
695:       ymax=-HUGE(ymax)
696:       xmint=HUGE(xmint)
697:       xmaxt=-HUGE(xmaxt)
698:       ymint=HUGE(ymint)
699:       ymaxt=-HUGE(ymaxt)
700: !-----------------------------------------------------------------------
701: !     get arrays and sizes.
702: !-----------------------------------------------------------------------
703:       CALL VecGetArray(v,va,iv,ierr)
704:       CALL VecGetSize(v,vn,ierr)
705:       CALL VecGetSize(fa%l,ln,ierr)
706: !-----------------------------------------------------------------------
707: !     fill arrays.
708: !-----------------------------------------------------------------------
709:       DO j=0,2
710:          IF (vn == ln)THEN
711:             offset=iv+2*fa%offl(j)
712:             patch(j)%mx=fa%ml(j)-1
713:             patch(j)%my=fa%nl(j)-1
714:          ELSE
715:             offset=iv+2*fa%offg(j)
716:             patch(j)%mx=fa%mg(j)-1
717:             patch(j)%my=fa%ng(j)-1
718:          ENDIF
719:          ALLOCATE(patch(j)%xy(2,0:patch(j)%mx,0:patch(j)%my))
720:          patch(j)%xy=RESHAPE(va(offset+1:offset+SIZE(patch(j)%xy)),                &
721:      &        (/2,patch(j)%mx+1,patch(j)%my+1/))
722:       ENDDO
723: !-----------------------------------------------------------------------
724: !     compute extrema over processor.
725: !-----------------------------------------------------------------------
726:       DO j=0,2
727:          xmin=MIN(xmin,MINVAL(patch(j)%xy(1,:,:)))
728:          xmax=MAX(xmax,MAXVAL(patch(j)%xy(1,:,:)))
729:          ymin=MIN(ymin,MINVAL(patch(j)%xy(2,:,:)))
730:          ymax=MAX(ymax,MAXVAL(patch(j)%xy(2,:,:)))
731:       ENDDO
732: !-----------------------------------------------------------------------
733: !     compute global extrema.
734: !-----------------------------------------------------------------------
735:       CALL MPI_Allreduce(xmin,xmint,1,MPI_DOUBLE_PRECISION,MPI_MIN,                 &
736:      &     PETSC_COMM_WORLD,ierr)
737:       CALL MPI_Allreduce(xmax,xmaxt,1,MPI_DOUBLE_PRECISION,MPI_MAX,                 &
738:      &     PETSC_COMM_WORLD,ierr)
739:       CALL MPI_Allreduce(ymin,ymint,1,MPI_DOUBLE_PRECISION,MPI_MIN,                 &
740:      &     PETSC_COMM_WORLD,ierr)
741:       CALL MPI_Allreduce(ymax,ymaxt,1,MPI_DOUBLE_PRECISION,MPI_MAX,                 &
742:      &     PETSC_COMM_WORLD,ierr)
743: !-----------------------------------------------------------------------
744: !     set margins.
745: !-----------------------------------------------------------------------
746:       xmin=xmint-.2*(xmaxt-xmint)
747:       xmax=xmaxt+.2*(xmaxt-xmint)
748:       ymin=ymint-.2*(ymaxt-ymint)
749:       ymax=ymaxt+.2*(ymaxt-ymint)
750: !-----------------------------------------------------------------------
751: !     draw it.
752: !-----------------------------------------------------------------------
753: #if defined(PETSC_HAVE_X) || defined(PETSC_HAVE_OPENGL)
754:       CALL PetscDrawCreate(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE,                  &
755:      &     PETSC_DECIDE,700,700,draw,ierr)
756:       CALL PetscDrawSetFromOptions(draw,ierr)
757:       CALL PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax,ierr)
758:       CALL PetscDrawZoom(draw,barry_draw_patch,patch,ierr)
759: #endif
760: !-----------------------------------------------------------------------
761: !     clean up.
762: !-----------------------------------------------------------------------
763:       CALL VecRestoreArray(v,va,iv,ierr)
764: #if defined(PETSC_HAVE_X)
765:       CALL PetscDrawDestroy(draw,ierr)
766: #endif
767: !-----------------------------------------------------------------------
768: !     terminate.
769: !-----------------------------------------------------------------------
770:       RETURN
771:       END SUBROUTINE barry_draw_fa
772: !-----------------------------------------------------------------------
773: !     subprogram 9. barry_map_region3.
774: !     fills region 3 coordinates.
775: !-----------------------------------------------------------------------
776: !-----------------------------------------------------------------------
777: !     declarations.
778: !-----------------------------------------------------------------------
779:       SUBROUTINE barry_map_region3(fa,g)

781:       TYPE(fa_type), INTENT(INOUT) :: fa
782:       Vec, INTENT(INOUT) :: g

784:       INTEGER :: i,j,k,x,y,m,n,ig
785:       REAL(8) :: r0=1,theta0=pi/6,r,theta,dr,dt
786:       REAL(8), DIMENSION(1) :: ga
787: !-----------------------------------------------------------------------
788: !     format statements.
789: !-----------------------------------------------------------------------
790:  10   FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
791:  20   FORMAT(2i6,4f11.3)
792: !-----------------------------------------------------------------------
793: !     preliminary computations.
794: !-----------------------------------------------------------------------
795:       dr=r0/(fa%r2-1)
796:       dt=twopi/(3*(fa%p1-fa%p2-1))
797:       CALL barry_get_global_corners(fa,2,x,y,m,n)
798: !-----------------------------------------------------------------------
799: !     fill array.
800: !-----------------------------------------------------------------------
801:       CALL VecGetArray(g,ga,ig,ierr)
802:       k=ig+2*fa%offg(2)
803:       IF (diagnose)THEN
804:          WRITE(*,'(a)')"region 3"
805:          WRITE(*,10)
806:       ENDIF
807:       DO j=y,y+n-1
808:          r=r0+j*dr
809:          DO i=x,x+m-1
810:             theta=theta0+i*dt
811:             ga(k+1)=r*COS(theta)
812:             ga(k+2)=r*SIN(theta)-4*r0
813:             IF (diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
814:             k=k+2
815:          ENDDO
816:          IF (diagnose)WRITE(*,10)
817:       ENDDO
818:       CALL VecRestoreArray(g,ga,ig,ierr)
819: !-----------------------------------------------------------------------
820: !     terminate.
821: !-----------------------------------------------------------------------
822:       RETURN
823:       END SUBROUTINE barry_map_region3
824: !-----------------------------------------------------------------------
825: !     subprogram 10. barry_map_region2.
826: !     fills region 2 coordinates.
827: !-----------------------------------------------------------------------
828: !-----------------------------------------------------------------------
829: !     declarations.
830: !-----------------------------------------------------------------------
831:       SUBROUTINE barry_map_region2(fa,g)

833:       TYPE(fa_type), INTENT(INOUT) :: fa
834:       Vec, INTENT(INOUT) :: g

836:       INTEGER :: i,j,k,x,y,m,n,ig
837:       REAL(8) :: r0=1,theta0=-pi/2,r,theta,dr,dt
838:       REAL(8), DIMENSION(1) :: ga
839: !-----------------------------------------------------------------------
840: !     format statements.
841: !-----------------------------------------------------------------------
842:  10   FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
843:  20   FORMAT(2i6,4f11.3)
844: !-----------------------------------------------------------------------
845: !     preliminary computations.
846: !-----------------------------------------------------------------------
847:       dr=r0/(fa%r2-1)
848:       dt=twopi/fa%p2
849:       CALL barry_get_global_corners(fa,1,x,y,m,n)
850: !-----------------------------------------------------------------------
851: !     fill array.
852: !-----------------------------------------------------------------------
853:       CALL VecGetArray(g,ga,ig,ierr)
854:       k=ig+2*fa%offg(1)
855:       IF (diagnose)THEN
856:          WRITE(*,'(a)')"region 2"
857:          WRITE(*,10)
858:       ENDIF
859:       DO j=y,y+n-1
860:          r=r0+j*dr
861:          DO i=x,x+m-1
862:             theta=theta0+i*dt
863:             ga(k+1)=r*COS(theta)
864:             ga(k+2)=r*SIN(theta)
865:             IF (diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
866:             k=k+2
867:          ENDDO
868:          IF (diagnose)WRITE(*,10)
869:       ENDDO
870:       CALL VecRestoreArray(g,ga,ig,ierr)
871: !-----------------------------------------------------------------------
872: !     terminate.
873: !-----------------------------------------------------------------------
874:       RETURN
875:       END SUBROUTINE barry_map_region2
876: !-----------------------------------------------------------------------
877: !     subprogram 11. barry_map_region1.
878: !     fills region 1 coordinates.
879: !-----------------------------------------------------------------------
880: !-----------------------------------------------------------------------
881: !     declarations.
882: !-----------------------------------------------------------------------
883:       SUBROUTINE barry_map_region1(fa,g)

885:       TYPE(fa_type), INTENT(INOUT) :: fa
886:       Vec, INTENT(INOUT) :: g

888:       INTEGER :: i,j,k,x,y,m,n,ig
889:       REAL(8) :: r0=1,r,theta,dr,dt1,dt3
890:       REAL(8), DIMENSION(1) :: ga
891: !-----------------------------------------------------------------------
892: !     format statements.
893: !-----------------------------------------------------------------------
894:  10   FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
895:  20   FORMAT(2i6,4f11.3)
896: !-----------------------------------------------------------------------
897: !     preliminary computations.
898: !-----------------------------------------------------------------------
899:       dr=r0/(fa%r1-1)
900:       dt1=twopi/fa%p2
901:       dt3=twopi/(3*(fa%p1-fa%p2-1))
902:       CALL barry_get_global_corners(fa,0,x,y,m,n)
903: !-----------------------------------------------------------------------
904: !     fill array.
905: !-----------------------------------------------------------------------
906:       CALL VecGetArray(g,ga,ig,ierr)
907:       k=ig+2*fa%offg(0)
908:       IF (diagnose)THEN
909:          WRITE(*,'(a)')"region 1"
910:          WRITE(*,10)
911:       ENDIF
912:       DO j=y,y+n-1
913:          r=2*r0+j*dr
914:          DO i=x,x+m-1
915:             IF (i < (fa%p1-fa%p2)/2)THEN
916:                theta=i*dt3
917:                ga(k+1)=r*COS(theta)
918:                ga(k+2)=r*SIN(theta)-4*r0
919:             ELSEIF (i > fa%p2 + (fa%p1-fa%p2)/2)then
920:                theta=pi+i*dt3
921:                ga(k+1)=r*COS(PETSC_PI+i*dt3)
922:                ga(k+2)=r*SIN(PETSC_PI+i*dt3)-4*r0
923:             ELSE
924:                theta=(i-(fa%p1-fa%p2)/2)*dt1-pi/2
925:                ga(k+1)=r*COS(theta)
926:                ga(k+2)=r*SIN(theta)
927:             ENDIF
928:             IF (diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
929:             k=k+2
930:          ENDDO
931:          IF (diagnose)WRITE(*,10)
932:       ENDDO
933:       CALL VecRestoreArray(g,ga,ig,ierr)
934: !-----------------------------------------------------------------------
935: !     terminate.
936: !-----------------------------------------------------------------------
937:       RETURN
938:       END SUBROUTINE barry_map_region1
939:       END MODULE barry_mod
940: !-----------------------------------------------------------------------
941: !     subprogram 12. barry_main.
942: !     main program.
943: !-----------------------------------------------------------------------
944: !-----------------------------------------------------------------------
945: !     declarations.
946: !-----------------------------------------------------------------------
947:       PROGRAM barry_main
948:       USE barry_mod
949:       IMPLICIT NONE

951:       TYPE(fa_type) :: fa
952:       Vec :: g,l
953: !-----------------------------------------------------------------------
954: !     initialize and create structure, and get vectors.
955: !-----------------------------------------------------------------------
956:       CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
957:       CALL barry_create_fa(fa)
958:       CALL barry_get_global_vector(fa,g)
959:       CALL barry_get_local_vector(fa,l)
960: !-----------------------------------------------------------------------
961: !     map regions.
962: !-----------------------------------------------------------------------
963:       CALL barry_map_region1(fa,g)
964:       CALL barry_map_region2(fa,g)
965:       CALL barry_map_region3(fa,g)
966: !-----------------------------------------------------------------------
967: !     graphic output.
968: !-----------------------------------------------------------------------
969:       CALL barry_global_to_local(fa,g,l)
970:       CALL barry_draw_fa(fa,g)
971:       CALL barry_draw_fa(fa,l)
972: !-----------------------------------------------------------------------
973: !     clean up and finalize.
974: !-----------------------------------------------------------------------
975:       CALL VecDestroy(g,ierr)
976:       CALL VecDestroy(l,ierr)
977:       CALL barry_destroy_fa(fa)
978:       CALL PetscFinalize(ierr)
979: !-----------------------------------------------------------------------
980: !     terminate.
981: !-----------------------------------------------------------------------
982:       END PROGRAM barry_main