Actual source code: ex6f90.F

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

628:       PetscDraw, INTENT(IN) :: draw
629:       TYPE(patch_type), DIMENSION(0:2), INTENT(IN) :: patch
630:       INTEGER, INTENT(OUT) :: ierr

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

674:       TYPE(fa_type), INTENT(IN) :: fa
675:       Vec, INTENT(IN) :: v

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

774:       TYPE(fa_type), INTENT(INOUT) :: fa
775:       Vec, INTENT(INOUT) :: g

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

826:       TYPE(fa_type), INTENT(INOUT) :: fa
827:       Vec, INTENT(INOUT) :: g

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

878:       TYPE(fa_type), INTENT(INOUT) :: fa
879:       Vec, INTENT(INOUT) :: g

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

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