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(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD, &
148: & ierr)
149: CALL VecScatterEnd(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD, &
150: & 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(vscat,globalvec,tovec,INSERT_VALUES, &
510: & SCATTER_FORWARD,ierr)
511: CALL VecScatterEnd(vscat,globalvec,tovec,INSERT_VALUES, &
512: & SCATTER_FORWARD,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: #if defined(PETSC_HAVE_X11)
749: CALL PetscDrawOpenX(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE, &
750: & PETSC_DECIDE,700,700,draw,ierr)
751: CALL PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax,ierr)
752: CALL PetscDrawZoom(draw,barry_draw_patch,patch,ierr)
753: #endif
754: !-----------------------------------------------------------------------
755: ! clean up.
756: !-----------------------------------------------------------------------
757: CALL VecRestoreArray(v,va,iv,ierr)
758: #if defined(PETSC_HAVE_X11)
759: CALL PetscDrawDestroy(draw,ierr)
760: #endif
761: !-----------------------------------------------------------------------
762: ! terminate.
763: !-----------------------------------------------------------------------
764: RETURN
765: END SUBROUTINE barry_draw_fa
766: !-----------------------------------------------------------------------
767: ! subprogram 9. barry_map_region3.
768: ! fills region 3 coordinates.
769: !-----------------------------------------------------------------------
770: !-----------------------------------------------------------------------
771: ! declarations.
772: !-----------------------------------------------------------------------
773: SUBROUTINE barry_map_region3(fa,g)
775: TYPE(fa_type), INTENT(INOUT) :: fa
776: Vec, INTENT(INOUT) :: g
778: INTEGER :: i,j,k,x,y,m,n,ig
779: REAL(8) :: r0=1,theta0=pi/6,r,theta,dr,dt
780: REAL(8), DIMENSION(1) :: ga
781: !-----------------------------------------------------------------------
782: ! format statements.
783: !-----------------------------------------------------------------------
784: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
785: 20 FORMAT(2i6,4f11.3)
786: !-----------------------------------------------------------------------
787: ! preliminary computations.
788: !-----------------------------------------------------------------------
789: dr=r0/(fa%r2-1)
790: dt=twopi/(3*(fa%p1-fa%p2-1))
791: CALL barry_get_global_corners(fa,2,x,y,m,n)
792: !-----------------------------------------------------------------------
793: ! fill array.
794: !-----------------------------------------------------------------------
795: CALL VecGetArray(g,ga,ig,ierr)
796: k=ig+2*fa%offg(2)
797: IF(diagnose)THEN
798: WRITE(*,'(a)')"region 3"
799: WRITE(*,10)
800: ENDIF
801: DO j=y,y+n-1
802: r=r0+j*dr
803: DO i=x,x+m-1
804: theta=theta0+i*dt
805: ga(k+1)=r*COS(theta)
806: ga(k+2)=r*SIN(theta)-4*r0
807: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
808: k=k+2
809: ENDDO
810: IF(diagnose)WRITE(*,10)
811: ENDDO
812: CALL VecRestoreArray(g,ga,ig,ierr)
813: !-----------------------------------------------------------------------
814: ! terminate.
815: !-----------------------------------------------------------------------
816: RETURN
817: END SUBROUTINE barry_map_region3
818: !-----------------------------------------------------------------------
819: ! subprogram 10. barry_map_region2.
820: ! fills region 2 coordinates.
821: !-----------------------------------------------------------------------
822: !-----------------------------------------------------------------------
823: ! declarations.
824: !-----------------------------------------------------------------------
825: SUBROUTINE barry_map_region2(fa,g)
827: TYPE(fa_type), INTENT(INOUT) :: fa
828: Vec, INTENT(INOUT) :: g
830: INTEGER :: i,j,k,x,y,m,n,ig
831: REAL(8) :: r0=1,theta0=-pi/2,r,theta,dr,dt
832: REAL(8), DIMENSION(1) :: ga
833: !-----------------------------------------------------------------------
834: ! format statements.
835: !-----------------------------------------------------------------------
836: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
837: 20 FORMAT(2i6,4f11.3)
838: !-----------------------------------------------------------------------
839: ! preliminary computations.
840: !-----------------------------------------------------------------------
841: dr=r0/(fa%r2-1)
842: dt=twopi/fa%p2
843: CALL barry_get_global_corners(fa,1,x,y,m,n)
844: !-----------------------------------------------------------------------
845: ! fill array.
846: !-----------------------------------------------------------------------
847: CALL VecGetArray(g,ga,ig,ierr)
848: k=ig+2*fa%offg(1)
849: IF(diagnose)THEN
850: WRITE(*,'(a)')"region 2"
851: WRITE(*,10)
852: ENDIF
853: DO j=y,y+n-1
854: r=r0+j*dr
855: DO i=x,x+m-1
856: theta=theta0+i*dt
857: ga(k+1)=r*COS(theta)
858: ga(k+2)=r*SIN(theta)
859: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
860: k=k+2
861: ENDDO
862: IF(diagnose)WRITE(*,10)
863: ENDDO
864: CALL VecRestoreArray(g,ga,ig,ierr)
865: !-----------------------------------------------------------------------
866: ! terminate.
867: !-----------------------------------------------------------------------
868: RETURN
869: END SUBROUTINE barry_map_region2
870: !-----------------------------------------------------------------------
871: ! subprogram 11. barry_map_region1.
872: ! fills region 1 coordinates.
873: !-----------------------------------------------------------------------
874: !-----------------------------------------------------------------------
875: ! declarations.
876: !-----------------------------------------------------------------------
877: SUBROUTINE barry_map_region1(fa,g)
879: TYPE(fa_type), INTENT(INOUT) :: fa
880: Vec, INTENT(INOUT) :: g
882: INTEGER :: i,j,k,x,y,m,n,ig
883: REAL(8) :: r0=1,r,theta,dr,dt1,dt3
884: REAL(8), DIMENSION(1) :: ga
885: !-----------------------------------------------------------------------
886: ! format statements.
887: !-----------------------------------------------------------------------
888: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
889: 20 FORMAT(2i6,4f11.3)
890: !-----------------------------------------------------------------------
891: ! preliminary computations.
892: !-----------------------------------------------------------------------
893: dr=r0/(fa%r1-1)
894: dt1=twopi/fa%p2
895: dt3=twopi/(3*(fa%p1-fa%p2-1))
896: CALL barry_get_global_corners(fa,0,x,y,m,n)
897: !-----------------------------------------------------------------------
898: ! fill array.
899: !-----------------------------------------------------------------------
900: CALL VecGetArray(g,ga,ig,ierr)
901: k=ig+2*fa%offg(0)
902: IF(diagnose)THEN
903: WRITE(*,'(a)')"region 1"
904: WRITE(*,10)
905: ENDIF
906: DO j=y,y+n-1
907: r=2*r0+j*dr
908: DO i=x,x+m-1
909: IF(i < (fa%p1-fa%p2)/2)THEN
910: theta=i*dt3
911: ga(k+1)=r*COS(theta)
912: ga(k+2)=r*SIN(theta)-4*r0
913: ELSEIF(i > fa%p2 + (fa%p1-fa%p2)/2)then
914: theta=pi+i*dt3
915: ga(k+1)=r*COS(PETSC_PI+i*dt3)
916: ga(k+2)=r*SIN(PETSC_PI+i*dt3)-4*r0
917: ELSE
918: theta=(i-(fa%p1-fa%p2)/2)*dt1-pi/2
919: ga(k+1)=r*COS(theta)
920: ga(k+2)=r*SIN(theta)
921: ENDIF
922: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
923: k=k+2
924: ENDDO
925: IF(diagnose)WRITE(*,10)
926: ENDDO
927: CALL VecRestoreArray(g,ga,ig,ierr)
928: !-----------------------------------------------------------------------
929: ! terminate.
930: !-----------------------------------------------------------------------
931: RETURN
932: END SUBROUTINE barry_map_region1
933: END MODULE barry_mod
934: !-----------------------------------------------------------------------
935: ! subprogram 12. barry_main.
936: ! main program.
937: !-----------------------------------------------------------------------
938: !-----------------------------------------------------------------------
939: ! declarations.
940: !-----------------------------------------------------------------------
941: PROGRAM barry_main
942: USE barry_mod
943: IMPLICIT NONE
945: TYPE(fa_type) :: fa
946: Vec :: g,l
947: !-----------------------------------------------------------------------
948: ! initialize and create structure, and get vectors.
949: !-----------------------------------------------------------------------
950: CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
951: CALL barry_create_fa(fa)
952: CALL barry_get_global_vector(fa,g)
953: CALL barry_get_local_vector(fa,l)
954: !-----------------------------------------------------------------------
955: ! map regions.
956: !-----------------------------------------------------------------------
957: CALL barry_map_region1(fa,g)
958: CALL barry_map_region2(fa,g)
959: CALL barry_map_region3(fa,g)
960: !-----------------------------------------------------------------------
961: ! graphic output.
962: !-----------------------------------------------------------------------
963: CALL barry_global_to_local(fa,g,l)
964: CALL barry_draw_fa(fa,g)
965: CALL barry_draw_fa(fa,l)
966: !-----------------------------------------------------------------------
967: ! clean up and finalize.
968: !-----------------------------------------------------------------------
969: CALL VecDestroy(g,ierr)
970: CALL VecDestroy(l,ierr)
971: CALL barry_destroy_fa(fa)
972: CALL PetscFinalize(ierr)
973: !-----------------------------------------------------------------------
974: ! terminate.
975: !-----------------------------------------------------------------------
976: END PROGRAM barry_main