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