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 include/finclude/petsc.h
 34:  #include include/finclude/petscsys.h
 35:  #include include/finclude/petscvec.h
 36:  #include include/finclude/petscda.h
 37:  #include include/finclude/petscmat.h
 38:  #include include/finclude/petscis.h
 39:  #include include/finclude/petscao.h
 40:  #include include/finclude/petscviewer.h
 41:  #include include/finclude/petscdraw.h

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

771:       TYPE(fa_type), INTENT(INOUT) :: fa
772:       Vec, INTENT(INOUT) :: g

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

823:       TYPE(fa_type), INTENT(INOUT) :: fa
824:       Vec, INTENT(INOUT) :: g

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

875:       TYPE(fa_type), INTENT(INOUT) :: fa
876:       Vec, INTENT(INOUT) :: g

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

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