Actual source code: ex6f90.F
petsc-3.3-p7 2013-05-11
1: !-----------------------------------------------------------------------
2: !
3: ! manages tokamak edge region.
4: !
5: ! This is a translation of ex6.c into F90 by Alan Glasser, LANL
6: !-----------------------------------------------------------------------
7: !-----------------------------------------------------------------------
8: ! code organization.
9: !-----------------------------------------------------------------------
10: ! 0. barry_mod.
11: ! 1. barry_get_global_corners.
12: ! 2. barry_get_global_vector.
13: ! 3. barry_get_local_vector.
14: ! 4. barry_global_to_local.
15: ! 5. barry_destroy_fa.
16: ! 6. barry_create_fa.
17: ! 7. barry_draw_patch.
18: ! 8. barry_draw_fa.
19: ! 9. barry_map_region3.
20: ! 10. barry_map_region2.
21: ! 11. barry_map_region1.
22: ! 12. barry_main.
23: !-----------------------------------------------------------------------
24: ! subprogram 0. barry_mod.
25: ! module declarations.
26: !-----------------------------------------------------------------------
27: !-----------------------------------------------------------------------
28: ! declarations.
29: !-----------------------------------------------------------------------
30: MODULE barry_mod
31: IMPLICIT NONE
33: #include <finclude/petscsys.h>
34: #include <finclude/petscvec.h>
35: #include <finclude/petscdmda.h>
36: #include <finclude/petscmat.h>
37: #include <finclude/petscis.h>
38: #include <finclude/petscao.h>
39: #include <finclude/petscviewer.h>
40: #include <finclude/petscdraw.h>
42: TYPE :: fa_type
43: MPI_Comm, DIMENSION(0:2) :: comm
44: INTEGER, DIMENSION(0:2) :: xl,yl,ml,nl,xg,yg,mg,ng,offl,offg
45: Vec :: g,l
46: VecScatter :: vscat
47: INTEGER :: p1,p2,r1,r2,r1g,r2g,sw
48: END TYPE fa_type
50: TYPE :: patch_type
51: INTEGER :: mx,my
52: PetscScalar, DIMENSION(:,:,:), POINTER :: xy
53: END TYPE patch_type
55: LOGICAL, PARAMETER :: diagnose=.FALSE.
56: INTEGER :: ierr
57: REAL(8), PARAMETER :: pi=3.1415926535897932385_8,twopi=2*pi
59: CONTAINS
60: !-----------------------------------------------------------------------
61: ! subprogram 1. barry_get_global_corners.
62: ! returns global corner data.
63: !-----------------------------------------------------------------------
64: !-----------------------------------------------------------------------
65: ! declarations.
66: !-----------------------------------------------------------------------
67: SUBROUTINE barry_get_global_corners(fa,j,x,y,m,n)
69: TYPE(fa_type), INTENT(IN) :: fa
70: INTEGER, INTENT(IN) :: j
71: INTEGER, INTENT(OUT) :: x,y,m,n
72: !-----------------------------------------------------------------------
73: ! computations.
74: !-----------------------------------------------------------------------
75: IF(fa%comm(j) /= 0)THEN
76: x=fa%xg(j)
77: y=fa%yg(j)
78: m=fa%mg(j)
79: n=fa%ng(j)
80: ELSE
81: x=0
82: y=0
83: m=0
84: n=0
85: ENDIF
86: !-----------------------------------------------------------------------
87: ! terminate.
88: !-----------------------------------------------------------------------
89: RETURN
90: END SUBROUTINE barry_get_global_corners
91: !-----------------------------------------------------------------------
92: ! subprogram 2. barry_get_global_vector.
93: ! returns local vector.
94: !-----------------------------------------------------------------------
95: !-----------------------------------------------------------------------
96: ! declarations.
97: !-----------------------------------------------------------------------
98: SUBROUTINE barry_get_global_vector(fa,v)
100: TYPE(fa_type), INTENT(IN) :: fa
101: Vec, INTENT(OUT) :: v
102: !-----------------------------------------------------------------------
103: ! computations.
104: !-----------------------------------------------------------------------
105: CALL VecDuplicate(fa%g,v,ierr)
106: !-----------------------------------------------------------------------
107: ! terminate.
108: !-----------------------------------------------------------------------
109: RETURN
110: END SUBROUTINE barry_get_global_vector
111: !-----------------------------------------------------------------------
112: ! subprogram 3. barry_get_local_vector.
113: ! returns local vector.
114: !-----------------------------------------------------------------------
115: !-----------------------------------------------------------------------
116: ! declarations.
117: !-----------------------------------------------------------------------
118: SUBROUTINE barry_get_local_vector(fa,v)
120: TYPE(fa_type), INTENT(IN) :: fa
121: Vec, INTENT(OUT) :: v
122: !-----------------------------------------------------------------------
123: ! computations.
124: !-----------------------------------------------------------------------
125: CALL VecDuplicate(fa%l,v,ierr)
126: !-----------------------------------------------------------------------
127: ! terminate.
128: !-----------------------------------------------------------------------
129: RETURN
130: END SUBROUTINE barry_get_local_vector
131: !-----------------------------------------------------------------------
132: ! subprogram 4. barry_global_to_local.
133: ! performs VecScatter.
134: !-----------------------------------------------------------------------
135: !-----------------------------------------------------------------------
136: ! declarations.
137: !-----------------------------------------------------------------------
138: SUBROUTINE barry_global_to_local(fa,g,l)
140: TYPE(fa_type), INTENT(IN) :: fa
141: Vec, INTENT(IN) :: g
142: Vec, INTENT(OUT) :: l
143: !-----------------------------------------------------------------------
144: ! computations.
145: !-----------------------------------------------------------------------
146: CALL VecScatterBegin(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD, &
147: & ierr)
148: CALL VecScatterEnd(fa%vscat,g,l,INSERT_VALUES,SCATTER_FORWARD, &
149: & ierr)
150: !-----------------------------------------------------------------------
151: ! terminate.
152: !-----------------------------------------------------------------------
153: RETURN
154: END SUBROUTINE barry_global_to_local
155: !-----------------------------------------------------------------------
156: ! subprogram 5. barry_destroy_fa.
157: ! destroys fa.
158: !-----------------------------------------------------------------------
159: !-----------------------------------------------------------------------
160: ! declarations.
161: !-----------------------------------------------------------------------
162: SUBROUTINE barry_destroy_fa(fa)
164: TYPE(fa_type), INTENT(OUT) :: fa
165: !-----------------------------------------------------------------------
166: ! computations.
167: !-----------------------------------------------------------------------
168: CALL VecDestroy(fa%g,ierr)
169: CALL VecDestroy(fa%l,ierr)
170: CALL VecScatterDestroy(fa%vscat,ierr)
171: !-----------------------------------------------------------------------
172: ! terminate.
173: !-----------------------------------------------------------------------
174: RETURN
175: END SUBROUTINE barry_destroy_fa
176: !-----------------------------------------------------------------------
177: ! subprogram 6. barry_create_fa.
178: ! creates fa.
179: !-----------------------------------------------------------------------
180: !-----------------------------------------------------------------------
181: ! declarations.
182: !-----------------------------------------------------------------------
183: SUBROUTINE barry_create_fa(fa)
185: TYPE(fa_type), INTENT(OUT) :: fa
187: INTEGER :: rank,tonglobal,globalrstart,x,nx,y,ny,i,j,k,ig, &
188: & fromnglobal,nscat,nlocal,cntl1,cntl2,cntl3,il,it
189: INTEGER, DIMENSION(0:2) :: offt
190: INTEGER, DIMENSION(:), POINTER :: tonatural,fromnatural, &
191: & to,from,indices
192: PetscScalar, DIMENSION(1) :: globalarray
193: PetscScalar, DIMENSION(1) :: localarray
194: PetscScalar, DIMENSION(1) :: toarray
196: DM :: da1=0,da2=0,da3=0
197: Vec :: vl1=0,vl2=0,vl3=0
198: Vec :: vg1=0,vg2=0,vg3=0
199: AO :: toao,globalao
200: IS :: tois,globalis,is
201: Vec :: tovec,globalvec,localvec
202: VecScatter :: vscat
203: !-----------------------------------------------------------------------
204: ! set integer members of fa.
205: !-----------------------------------------------------------------------
206: fa%p1=24
207: fa%p2=15
208: fa%r1=6
209: fa%r2=6
210: fa%sw=1
211: fa%r1g=fa%r1+fa%sw
212: fa%r2g=fa%r2+fa%sw
213: !-----------------------------------------------------------------------
214: ! set up communicators.
215: !-----------------------------------------------------------------------
216: CALL MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
217: fa%comm=PETSC_COMM_WORLD
218: !-----------------------------------------------------------------------
219: ! set up distributed arrays.
220: !-----------------------------------------------------------------------
221: IF(fa%comm(1) /= 0)THEN
222: CALL DMDACreate2d(fa%comm(1),DMDA_BOUNDARY_PERIODIC, &
223: & DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
224: & fa%p2,fa%r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw, &
225: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da2,ierr)
226: CALL DMGetLocalVector(da2,vl2,ierr)
227: CALL DMGetGlobalVector(da2,vg2,ierr)
228: ENDIF
229: IF(fa%comm(2) /= 0)THEN
230: CALL DMDACreate2d(fa%comm(2),DMDA_BOUNDARY_NONE, &
231: & DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
232: & fa%p1-fa%p2,fa%r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw, &
233: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da3,ierr)
234: CALL DMGetLocalVector(da3,vl3,ierr)
235: CALL DMGetGlobalVector(da3,vg3,ierr)
236: ENDIF
237: IF(fa%comm(0) /= 0)THEN
238: CALL DMDACreate2d(fa%comm(0),DMDA_BOUNDARY_NONE, &
239: & DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
240: & fa%p1,fa%r1g,PETSC_DECIDE,PETSC_DECIDE,1,fa%sw, &
241: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da1,ierr)
242: CALL DMGetLocalVector(da1,vl1,ierr)
243: CALL DMGetGlobalVector(da1,vg1,ierr)
244: ENDIF
245: !-----------------------------------------------------------------------
246: ! count the number of unknowns owned on each processor and determine
247: ! the starting point of each processors ownership
248: ! for global vector with redundancy.
249: !-----------------------------------------------------------------------
250: tonglobal = 0
251: IF(fa%comm(1) /= 0)THEN
252: CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
253: tonglobal=tonglobal+nx*ny
254: ENDIF
255: IF(fa%comm(2) /= 0)THEN
256: CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
257: tonglobal=tonglobal+nx*ny
258: ENDIF
259: IF(fa%comm(0) /= 0)THEN
260: CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
261: tonglobal=tonglobal+nx*ny
262: ENDIF
263: WRITE(*,'("[",i1,"]",a,i3)') &
264: & rank," Number of unknowns owned ",tonglobal
265: !-----------------------------------------------------------------------
266: ! get tonatural number for each node
267: !-----------------------------------------------------------------------
268: ALLOCATE(tonatural(0:tonglobal))
269: tonglobal=0
270: IF(fa%comm(1) /= 0)THEN
271: CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
272: DO j=0,ny-1
273: DO i=0,nx-1
274: tonatural(tonglobal)=(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
275: tonglobal=tonglobal+1
276: ENDDO
277: ENDDO
278: ENDIF
279: IF(fa%comm(2) /= 0)THEN
280: CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
281: DO j=0,ny-1
282: DO i=0,nx-1
283: IF(x+i < (fa%p1-fa%p2)/2)THEN
284: tonatural(tonglobal)=x+i+fa%p1*(y+j)
285: ELSE
286: tonatural(tonglobal)=fa%p2+x+i+fa%p1*(y+j)
287: ENDIF
288: tonglobal=tonglobal+1
289: ENDDO
290: ENDDO
291: ENDIF
292: IF(fa%comm(0) /= 0)THEN
293: CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
294: DO j=0,ny-1
295: DO i=0,nx-1
296: tonatural(tonglobal)=fa%p1*fa%r2g+x+i+fa%p1*(y+j)
297: tonglobal=tonglobal+1
298: ENDDO
299: ENDDO
300: ENDIF
301: !-----------------------------------------------------------------------
302: ! diagnose tonatural.
303: !-----------------------------------------------------------------------
304: IF(diagnose)THEN
305: WRITE(*,'(a,i3,a)')"tonglobal = ",tonglobal,", tonatural = "
306: WRITE(*,'(2i4)')(i,tonatural(i),i=0,tonglobal-1)
307: ENDIF
308: !-----------------------------------------------------------------------
309: ! create application ordering and deallocate tonatural.
310: !-----------------------------------------------------------------------
311: CALL AOCreateBasic(PETSC_COMM_WORLD,tonglobal,tonatural, &
312: & PETSC_NULL_INTEGER,toao,ierr)
313: DEALLOCATE(tonatural)
314: !-----------------------------------------------------------------------
315: ! count the number of unknowns owned on each processor and determine
316: ! the starting point of each processors ownership for global vector
317: ! without redundancy.
318: !-----------------------------------------------------------------------
319: fromnglobal=0
320: fa%offg(1)=0
321: offt(1)=0
322: IF(fa%comm(1) /= 0)THEN
323: CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
324: offt(2)=nx*ny
325: IF(y+ny == fa%r2g)ny=ny-1
326: fromnglobal=fromnglobal+nx*ny
327: fa%offg(2)=fromnglobal
328: ELSE
329: offt(2)=0
330: fa%offg(2)=0
331: ENDIF
332: IF(fa%comm(2) /= 0)THEN
333: CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
334: offt(0)=offt(2)+nx*ny
335: IF(y+ny == fa%r2g)ny=ny-1
336: fromnglobal=fromnglobal+nx*ny
337: fa%offg(0)=fromnglobal
338: ELSE
339: offt(0)=offt(2)
340: fa%offg(0)=fromnglobal
341: ENDIF
342: IF(fa%comm(0) /= 0)THEN
343: CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
344: IF(y == 0)ny=ny-1
345: fromnglobal=fromnglobal+nx*ny
346: ENDIF
347: CALL MPI_Scan(fromnglobal,globalrstart,1,MPI_INTEGER, &
348: & MPI_SUM,PETSC_COMM_WORLD,ierr)
349: globalrstart=globalrstart-fromnglobal
350: WRITE(*,'("[",i1,"]",a,i3)') &
351: & rank," Number of unknowns owned ",fromnglobal
352: !-----------------------------------------------------------------------
353: ! get fromnatural number for each node.
354: !-----------------------------------------------------------------------
355: ALLOCATE(fromnatural(0:fromnglobal))
356: fromnglobal=0
357: IF(fa%comm(1) /= 0)THEN
358: CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
359: IF(y+ny == fa%r2g)ny=ny-1
360: fa%xg(1)=x
361: fa%yg(1)=y
362: fa%mg(1)=nx
363: fa%ng(1)=ny
364: CALL DMDAGetGhostCorners(da2,fa%xl(1),fa%yl(1),0,fa%ml(1), &
365: & fa%nl(1),0,ierr)
366: DO j=0,ny-1
367: DO i=0,nx-1
368: fromnatural(fromnglobal) &
369: & =(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
370: fromnglobal=fromnglobal+1
371: ENDDO
372: ENDDO
373: ENDIF
374: IF(fa%comm(2) /= 0)THEN
375: CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
376: IF(y+ny == fa%r2g)ny=ny-1
377: fa%xg(2)=x
378: fa%yg(2)=y
379: fa%mg(2)=nx
380: fa%ng(2)=ny
381: CALL DMDAGetGhostCorners(da3,fa%xl(2),fa%yl(2),0,fa%ml(2), &
382: & fa%nl(2),0,ierr)
383: DO j=0,ny-1
384: DO i=0,nx-1
385: IF(x+i < (fa%p1-fa%p2)/2)THEN
386: fromnatural(fromnglobal)=x+i+fa%p1*(y+j)
387: ELSE
388: fromnatural(fromnglobal)=fa%p2+x+i+fa%p1*(y+j)
389: ENDIF
390: fromnglobal=fromnglobal+1
391: ENDDO
392: ENDDO
393: ENDIF
394: IF(fa%comm(0) /= 0)THEN
395: CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
396: IF(y == 0)THEN
397: ny=ny-1
398: ELSE
399: y=y-1
400: ENDIF
401: fa%xg(0)=x
402: fa%yg(0)=y
403: fa%mg(0)=nx
404: fa%ng(0)=ny
405: CALL DMDAGetGhostCorners(da1,fa%xl(0),fa%yl(0),0,fa%ml(0), &
406: & fa%nl(0),0,ierr)
407: DO j=0,ny-1
408: DO i=0,nx-1
409: fromnatural(fromnglobal)=fa%p1*fa%r2+x+i+fa%p1*(y+j)
410: fromnglobal=fromnglobal+1
411: ENDDO
412: ENDDO
413: ENDIF
414: !-----------------------------------------------------------------------
415: ! diagnose fromnatural.
416: !-----------------------------------------------------------------------
417: IF(diagnose)THEN
418: WRITE(*,'(a,i3,a)')"fromnglobal = ",fromnglobal &
419: & ,", fromnatural = "
420: WRITE(*,'(2i4)')(i,fromnatural(i),i=0,fromnglobal-1)
421: ENDIF
422: !-----------------------------------------------------------------------
423: ! create application ordering and deallocate fromnatural.
424: !-----------------------------------------------------------------------
425: CALL AOCreateBasic(PETSC_COMM_WORLD,fromnglobal,fromnatural, &
426: & PETSC_NULL_INTEGER,globalao,ierr)
427: DEALLOCATE(fromnatural)
428: !-----------------------------------------------------------------------
429: ! set up scatter that updates 1 from 2 and 3 and 3 and 2 from 1
430: !-----------------------------------------------------------------------
431: ALLOCATE(to(0:tonglobal),from(0:tonglobal))
432: nscat=0
433: IF(fa%comm(1) /= 0)THEN
434: CALL DMDAGetCorners(da2,x,y,0,nx,ny,0,ierr)
435: DO j=0,ny-1
436: DO i=0,nx-1
437: to(nscat)=(fa%p1-fa%p2)/2+x+i+fa%p1*(y+j)
438: from(nscat)=to(nscat)
439: nscat=nscat+1
440: ENDDO
441: ENDDO
442: ENDIF
443: IF(fa%comm(2) /= 0)THEN
444: CALL DMDAGetCorners(da3,x,y,0,nx,ny,0,ierr)
445: DO j=0,ny-1
446: DO i=0,nx-1
447: IF(x+i < (fa%p1-fa%p2)/2)THEN
448: to(nscat)=x+i+fa%p1*(y+j)
449: ELSE
450: to(nscat)=fa%p2+x+i+fa%p1*(y+j)
451: ENDIF
452: from(nscat)=to(nscat)
453: nscat=nscat+1
454: ENDDO
455: ENDDO
456: ENDIF
457: IF(fa%comm(0) /= 0)THEN
458: CALL DMDAGetCorners(da1,x,y,0,nx,ny,0,ierr)
459: DO j=0,ny-1
460: DO i=0,nx-1
461: to(nscat)=fa%p1*fa%r2g+x+i+fa%p1*(y+j)
462: from(nscat)=fa%p1*(fa%r2-1)+x+i+fa%p1*(y+j)
463: nscat=nscat+1
464: ENDDO
465: ENDDO
466: ENDIF
467: !-----------------------------------------------------------------------
468: ! diagnose to and from.
469: !-----------------------------------------------------------------------
470: IF(diagnose)THEN
471: WRITE(*,'(a,i3,a)')"nscat = ",nscat,", to, from = "
472: WRITE(*,'(3i4)')(i,to(i),from(i),i=0,nscat-1)
473: ENDIF
474: !-----------------------------------------------------------------------
475: ! create vecscatter.
476: !-----------------------------------------------------------------------
477: CALL AOApplicationToPetsc(toao,nscat,to,ierr)
478: CALL AOApplicationToPetsc(globalao,nscat,from,ierr)
479: CALL ISCreateGeneral(PETSC_COMM_WORLD,nscat,to,PETSC_COPY_VALUES, &
480: & tois,ierr)
481: CALL ISCreateGeneral(PETSC_COMM_WORLD,nscat,from,PETSC_COPY_VALUES &
482: & ,globalis,ierr)
483: CALL VecCreateMPI(PETSC_COMM_WORLD,tonglobal,PETSC_DETERMINE, &
484: & tovec,ierr)
485: CALL VecCreateMPI(PETSC_COMM_WORLD,fromnglobal,PETSC_DETERMINE, &
486: & globalvec,ierr)
487: CALL VecScatterCreate(globalvec,globalis,tovec,tois,vscat,ierr)
488: !-----------------------------------------------------------------------
489: ! clean up.
490: !-----------------------------------------------------------------------
491: CALL ISDestroy(tois,ierr)
492: CALL ISDestroy(globalis,ierr)
493: CALL AODestroy(globalao,ierr)
494: CALL AODestroy(toao,ierr)
495: DEALLOCATE(to,from)
496: !-----------------------------------------------------------------------
497: ! fill up global vector without redundant values with PETSc global
498: ! numbering
499: !-----------------------------------------------------------------------
500: CALL VecGetArray(globalvec,globalarray,ig,ierr)
501: k=ig
502: IF(diagnose)WRITE(*,'(a,i3,a)') &
503: & "fromnglobal = ",fromnglobal,", globalarray = "
504: DO i=0,fromnglobal-1
505: k=k+1
506: globalarray(k)=globalrstart+i
507: IF(diagnose)WRITE(*,'(i4,f11.3)')i,globalarray(k)
508: ENDDO
509: CALL VecRestoreArray(globalvec,globalarray,ig,ierr)
510: !-----------------------------------------------------------------------
511: ! scatter PETSc global indices to redundant valued array.
512: !-----------------------------------------------------------------------
513: CALL VecScatterBegin(vscat,globalvec,tovec,INSERT_VALUES, &
514: & SCATTER_FORWARD,ierr)
515: CALL VecScatterEnd(vscat,globalvec,tovec,INSERT_VALUES, &
516: & SCATTER_FORWARD,ierr)
517: !-----------------------------------------------------------------------
518: ! create local vector as concatenation of local vectors
519: !-----------------------------------------------------------------------
520: nlocal=0
521: cntl1=0
522: cntl2=0
523: cntl3=0
524: IF(fa%comm(1) /= 0)THEN
525: CALL VecGetSize(vl2,cntl2,ierr)
526: nlocal=nlocal+cntl2
527: ENDIF
528: IF(fa%comm(2) /= 0)THEN
529: CALL VecGetSize(vl3,cntl3,ierr)
530: nlocal=nlocal+cntl3
531: ENDIF
532: IF(fa%comm(0) /= 0)THEN
533: CALL VecGetSize(vl1,cntl1,ierr)
534: nlocal=nlocal+cntl1
535: ENDIF
536: fa%offl(0)=cntl2+cntl3
537: fa%offl(1)=0
538: fa%offl(2)=cntl2
539: CALL VecCreateSeq(PETSC_COMM_SELF,nlocal,localvec,ierr)
540: !-----------------------------------------------------------------------
541: ! cheat so that vl1, vl2, vl3 shared array memory with localvec.
542: !-----------------------------------------------------------------------
543: CALL VecGetArray(localvec,localarray,il,ierr)
544: CALL VecGetArray(tovec,toarray,it,ierr)
545: IF(fa%comm(1) /= 0)THEN
546: CALL VecPlaceArray(vl2,localarray(il+1+fa%offl(1)),ierr)
547: CALL VecPlaceArray(vg2,toarray(it+1+offt(1)),ierr)
548: CALL DMGlobalToLocalBegin(da2,vg2,INSERT_VALUES,vl2,ierr)
549: CALL DMGlobalToLocalEnd(da2,vg2,INSERT_VALUES,vl2,ierr)
550: CALL DMRestoreGlobalVector(da2,vg2,ierr)
551: ENDIF
552: IF(fa%comm(2) /= 0)THEN
553: CALL VecPlaceArray(vl3,localarray(il+1+fa%offl(2)),ierr)
554: CALL VecPlaceArray(vg3,toarray(it+1+offt(2)),ierr)
555: CALL DMGlobalToLocalBegin(da3,vg3,INSERT_VALUES,vl3,ierr)
556: CALL DMGlobalToLocalEnd(da3,vg3,INSERT_VALUES,vl3,ierr)
557: CALL DMRestoreGlobalVector(da3,vg3,ierr)
558: ENDIF
559: IF(fa%comm(0) /= 0)THEN
560: CALL VecPlaceArray(vl1,localarray(il+1+fa%offl(0)),ierr)
561: CALL VecPlaceArray(vg1,toarray(it+1+offt(0)),ierr)
562: CALL DMGlobalToLocalBegin(da1,vg1,INSERT_VALUES,vl1,ierr)
563: CALL DMGlobalToLocalEnd(da1,vg1,INSERT_VALUES,vl1,ierr)
564: CALL DMRestoreGlobalVector(da1,vg1,ierr)
565: ENDIF
566: CALL VecRestoreArray(localvec,localarray,il,ierr)
567: CALL VecRestoreArray(tovec,toarray,it,ierr)
568: !-----------------------------------------------------------------------
569: ! clean up.
570: !-----------------------------------------------------------------------
571: CALL VecScatterDestroy(vscat,ierr)
572: CALL VecDestroy(tovec,ierr)
573: !-----------------------------------------------------------------------
574: ! compute index set.
575: !-----------------------------------------------------------------------
576: ALLOCATE(indices(0:nlocal-1))
577: CALL VecGetArray(localvec,localarray,il,ierr)
578: k=il
579: IF(diagnose)WRITE(*,'(a,i3,a3,i4,a)') &
580: & "nlocal = ", nlocal,", offl = ",fa%offl,", indices = "
581: DO i=0,nlocal-1
582: k=k+1
583: indices(i)= int(localarray(k))
584: IF(diagnose)WRITE(*,'(2i4)')i,indices(i)
585: ENDDO
586: CALL VecRestoreArray(localvec,localarray,il,ierr)
587: CALL ISCreateBlock(PETSC_COMM_WORLD,2,nlocal,indices, &
588: & PETSC_COPY_VALUES,is,ierr)
589: DEALLOCATE(indices)
590: !-----------------------------------------------------------------------
591: ! create local and global vectors.
592: !-----------------------------------------------------------------------
593: CALL VecCreateSeq(PETSC_COMM_SELF,2*nlocal,fa%l,ierr)
594: CALL VecCreateMPI(PETSC_COMM_WORLD,2*fromnglobal,PETSC_DETERMINE, &
595: & fa%g,ierr)
596: !-----------------------------------------------------------------------
597: ! create final scatter that goes directly from globalvec to localvec
598: ! this is the one to be used in the application code.
599: !-----------------------------------------------------------------------
600: CALL VecScatterCreate(fa%g,is,fa%l,PETSC_NULL_OBJECT, &
601: & fa%vscat,ierr)
602: !-----------------------------------------------------------------------
603: ! clean up.
604: !-----------------------------------------------------------------------
605: CALL ISDestroy(is,ierr)
606: CALL VecDestroy(globalvec,ierr)
607: CALL VecDestroy(localvec,ierr)
608: IF(fa%comm(0) /= 0)THEN
609: CALL DMRestoreLocalVector(da1,vl1,ierr)
610: CALL DMDestroy(da1,ierr)
611: ENDIF
612: IF(fa%comm(1) /= 0)THEN
613: CALL DMRestoreLocalVector(da2,vl2,ierr)
614: CALL DMDestroy(da2,ierr)
615: ENDIF
616: IF(fa%comm(2) /= 0)THEN
617: CALL DMRestoreLocalVector(da3,vl3,ierr)
618: CALL DMDestroy(da3,ierr)
619: ENDIF
620: !-----------------------------------------------------------------------
621: ! terminate.
622: !-----------------------------------------------------------------------
623: RETURN
624: END SUBROUTINE barry_create_fa
625: !-----------------------------------------------------------------------
626: ! subprogram 7. barry_draw_patch.
627: ! crude graphics to test that the ghost points are properly updated.
628: !-----------------------------------------------------------------------
629: !-----------------------------------------------------------------------
630: ! declarations.
631: !-----------------------------------------------------------------------
632: SUBROUTINE barry_draw_patch(draw,patch,ierr)
634: PetscDraw, INTENT(IN) :: draw
635: TYPE(patch_type), DIMENSION(0:2), INTENT(IN) :: patch
636: INTEGER, INTENT(OUT) :: ierr
638: INTEGER :: ix,iy,j
639: PetscReal, DIMENSION(4) :: x,y
640: !-----------------------------------------------------------------------
641: ! draw it.
642: !-----------------------------------------------------------------------
643: DO j=0,2
644: DO iy=1,patch(j)%my
645: DO ix=1,patch(j)%mx
646: x(1)=patch(j)%xy(1,ix-1,iy-1)
647: y(1)=patch(j)%xy(2,ix-1,iy-1)
648: x(2)=patch(j)%xy(1,ix,iy-1)
649: y(2)=patch(j)%xy(2,ix,iy-1)
650: x(3)=patch(j)%xy(1,ix,iy)
651: y(3)=patch(j)%xy(2,ix,iy)
652: x(4)=patch(j)%xy(1,ix-1,iy)
653: y(4)=patch(j)%xy(2,ix-1,iy)
654: CALL PetscDrawLine(draw,x(1),y(1),x(2),y(2), &
655: & PETSC_DRAW_BLACK,ierr)
656: CALL PetscDrawLine(draw,x(2),y(2),x(3),y(3), &
657: & PETSC_DRAW_BLACK,ierr)
658: CALL PetscDrawLine(draw,x(3),y(3),x(4),y(4), &
659: & PETSC_DRAW_BLACK,ierr)
660: CALL PetscDrawLine(draw,x(4),y(4),x(1),y(1), &
661: & PETSC_DRAW_BLACK,ierr)
662: ENDDO
663: ENDDO
664: ENDDO
665: !-----------------------------------------------------------------------
666: ! terminate.
667: !-----------------------------------------------------------------------
668: ierr=0
669: RETURN
670: END SUBROUTINE barry_draw_patch
671: !-----------------------------------------------------------------------
672: ! subprogram 8. barry_draw_fa.
673: ! deallocates local array.
674: !-----------------------------------------------------------------------
675: !-----------------------------------------------------------------------
676: ! declarations.
677: !-----------------------------------------------------------------------
678: SUBROUTINE barry_draw_fa(fa,v)
680: TYPE(fa_type), INTENT(IN) :: fa
681: Vec, INTENT(IN) :: v
683: INTEGER :: iv,vn,ln,j,offset
684: REAL(8), DIMENSION(1) :: va
685: TYPE(patch_type), DIMENSION(0:2) :: patch
686: PetscDraw :: draw
687: PetscReal :: xmin,xmax,ymin,ymax
688: PetscReal :: ymint,xmint,ymaxt,xmaxt
689: !-----------------------------------------------------------------------
690: ! set extrema.
691: !-----------------------------------------------------------------------
692: xmin=HUGE(xmin)
693: xmax=-HUGE(xmax)
694: ymin=HUGE(ymin)
695: ymax=-HUGE(ymax)
696: xmint=HUGE(xmint)
697: xmaxt=-HUGE(xmaxt)
698: ymint=HUGE(ymint)
699: ymaxt=-HUGE(ymaxt)
700: !-----------------------------------------------------------------------
701: ! get arrays and sizes.
702: !-----------------------------------------------------------------------
703: CALL VecGetArray(v,va,iv,ierr)
704: CALL VecGetSize(v,vn,ierr)
705: CALL VecGetSize(fa%l,ln,ierr)
706: !-----------------------------------------------------------------------
707: ! fill arrays.
708: !-----------------------------------------------------------------------
709: DO j=0,2
710: IF(vn == ln)THEN
711: offset=iv+2*fa%offl(j)
712: patch(j)%mx=fa%ml(j)-1
713: patch(j)%my=fa%nl(j)-1
714: ELSE
715: offset=iv+2*fa%offg(j)
716: patch(j)%mx=fa%mg(j)-1
717: patch(j)%my=fa%ng(j)-1
718: ENDIF
719: ALLOCATE(patch(j)%xy(2,0:patch(j)%mx,0:patch(j)%my))
720: patch(j)%xy=RESHAPE(va(offset+1:offset+SIZE(patch(j)%xy)), &
721: & (/2,patch(j)%mx+1,patch(j)%my+1/))
722: ENDDO
723: !-----------------------------------------------------------------------
724: ! compute extrema over processor.
725: !-----------------------------------------------------------------------
726: DO j=0,2
727: xmin=MIN(xmin,MINVAL(patch(j)%xy(1,:,:)))
728: xmax=MAX(xmax,MAXVAL(patch(j)%xy(1,:,:)))
729: ymin=MIN(ymin,MINVAL(patch(j)%xy(2,:,:)))
730: ymax=MAX(ymax,MAXVAL(patch(j)%xy(2,:,:)))
731: ENDDO
732: !-----------------------------------------------------------------------
733: ! compute global extrema.
734: !-----------------------------------------------------------------------
735: CALL MPI_Allreduce(xmin,xmint,1,MPI_DOUBLE_PRECISION,MPI_MIN, &
736: & PETSC_COMM_WORLD,ierr)
737: CALL MPI_Allreduce(xmax,xmaxt,1,MPI_DOUBLE_PRECISION,MPI_MAX, &
738: & PETSC_COMM_WORLD,ierr)
739: CALL MPI_Allreduce(ymin,ymint,1,MPI_DOUBLE_PRECISION,MPI_MIN, &
740: & PETSC_COMM_WORLD,ierr)
741: CALL MPI_Allreduce(ymax,ymaxt,1,MPI_DOUBLE_PRECISION,MPI_MAX, &
742: & PETSC_COMM_WORLD,ierr)
743: !-----------------------------------------------------------------------
744: ! set margins.
745: !-----------------------------------------------------------------------
746: xmin=xmint-.2*(xmaxt-xmint)
747: xmax=xmaxt+.2*(xmaxt-xmint)
748: ymin=ymint-.2*(ymaxt-ymint)
749: ymax=ymaxt+.2*(ymaxt-ymint)
750: !-----------------------------------------------------------------------
751: ! draw it.
752: !-----------------------------------------------------------------------
753: #if defined(PETSC_HAVE_X)
754: CALL PetscDrawOpenX(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE, &
755: & PETSC_DECIDE,700,700,draw,ierr)
756: CALL PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax,ierr)
757: CALL PetscDrawZoom(draw,barry_draw_patch,patch,ierr)
758: #endif
759: !-----------------------------------------------------------------------
760: ! clean up.
761: !-----------------------------------------------------------------------
762: CALL VecRestoreArray(v,va,iv,ierr)
763: #if defined(PETSC_HAVE_X)
764: CALL PetscDrawDestroy(draw,ierr)
765: #endif
766: !-----------------------------------------------------------------------
767: ! terminate.
768: !-----------------------------------------------------------------------
769: RETURN
770: END SUBROUTINE barry_draw_fa
771: !-----------------------------------------------------------------------
772: ! subprogram 9. barry_map_region3.
773: ! fills region 3 coordinates.
774: !-----------------------------------------------------------------------
775: !-----------------------------------------------------------------------
776: ! declarations.
777: !-----------------------------------------------------------------------
778: SUBROUTINE barry_map_region3(fa,g)
780: TYPE(fa_type), INTENT(INOUT) :: fa
781: Vec, INTENT(INOUT) :: g
783: INTEGER :: i,j,k,x,y,m,n,ig
784: REAL(8) :: r0=1,theta0=pi/6,r,theta,dr,dt
785: REAL(8), DIMENSION(1) :: ga
786: !-----------------------------------------------------------------------
787: ! format statements.
788: !-----------------------------------------------------------------------
789: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
790: 20 FORMAT(2i6,4f11.3)
791: !-----------------------------------------------------------------------
792: ! preliminary computations.
793: !-----------------------------------------------------------------------
794: dr=r0/(fa%r2-1)
795: dt=twopi/(3*(fa%p1-fa%p2-1))
796: CALL barry_get_global_corners(fa,2,x,y,m,n)
797: !-----------------------------------------------------------------------
798: ! fill array.
799: !-----------------------------------------------------------------------
800: CALL VecGetArray(g,ga,ig,ierr)
801: k=ig+2*fa%offg(2)
802: IF(diagnose)THEN
803: WRITE(*,'(a)')"region 3"
804: WRITE(*,10)
805: ENDIF
806: DO j=y,y+n-1
807: r=r0+j*dr
808: DO i=x,x+m-1
809: theta=theta0+i*dt
810: ga(k+1)=r*COS(theta)
811: ga(k+2)=r*SIN(theta)-4*r0
812: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
813: k=k+2
814: ENDDO
815: IF(diagnose)WRITE(*,10)
816: ENDDO
817: CALL VecRestoreArray(g,ga,ig,ierr)
818: !-----------------------------------------------------------------------
819: ! terminate.
820: !-----------------------------------------------------------------------
821: RETURN
822: END SUBROUTINE barry_map_region3
823: !-----------------------------------------------------------------------
824: ! subprogram 10. barry_map_region2.
825: ! fills region 2 coordinates.
826: !-----------------------------------------------------------------------
827: !-----------------------------------------------------------------------
828: ! declarations.
829: !-----------------------------------------------------------------------
830: SUBROUTINE barry_map_region2(fa,g)
832: TYPE(fa_type), INTENT(INOUT) :: fa
833: Vec, INTENT(INOUT) :: g
835: INTEGER :: i,j,k,x,y,m,n,ig
836: REAL(8) :: r0=1,theta0=-pi/2,r,theta,dr,dt
837: REAL(8), DIMENSION(1) :: ga
838: !-----------------------------------------------------------------------
839: ! format statements.
840: !-----------------------------------------------------------------------
841: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
842: 20 FORMAT(2i6,4f11.3)
843: !-----------------------------------------------------------------------
844: ! preliminary computations.
845: !-----------------------------------------------------------------------
846: dr=r0/(fa%r2-1)
847: dt=twopi/fa%p2
848: CALL barry_get_global_corners(fa,1,x,y,m,n)
849: !-----------------------------------------------------------------------
850: ! fill array.
851: !-----------------------------------------------------------------------
852: CALL VecGetArray(g,ga,ig,ierr)
853: k=ig+2*fa%offg(1)
854: IF(diagnose)THEN
855: WRITE(*,'(a)')"region 2"
856: WRITE(*,10)
857: ENDIF
858: DO j=y,y+n-1
859: r=r0+j*dr
860: DO i=x,x+m-1
861: theta=theta0+i*dt
862: ga(k+1)=r*COS(theta)
863: ga(k+2)=r*SIN(theta)
864: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
865: k=k+2
866: ENDDO
867: IF(diagnose)WRITE(*,10)
868: ENDDO
869: CALL VecRestoreArray(g,ga,ig,ierr)
870: !-----------------------------------------------------------------------
871: ! terminate.
872: !-----------------------------------------------------------------------
873: RETURN
874: END SUBROUTINE barry_map_region2
875: !-----------------------------------------------------------------------
876: ! subprogram 11. barry_map_region1.
877: ! fills region 1 coordinates.
878: !-----------------------------------------------------------------------
879: !-----------------------------------------------------------------------
880: ! declarations.
881: !-----------------------------------------------------------------------
882: SUBROUTINE barry_map_region1(fa,g)
884: TYPE(fa_type), INTENT(INOUT) :: fa
885: Vec, INTENT(INOUT) :: g
887: INTEGER :: i,j,k,x,y,m,n,ig
888: REAL(8) :: r0=1,r,theta,dr,dt1,dt3
889: REAL(8), DIMENSION(1) :: ga
890: !-----------------------------------------------------------------------
891: ! format statements.
892: !-----------------------------------------------------------------------
893: 10 FORMAT(/5x,"j",5x,"i",6x,"r",8x,"theta",8x,"x",10x,"y"/)
894: 20 FORMAT(2i6,4f11.3)
895: !-----------------------------------------------------------------------
896: ! preliminary computations.
897: !-----------------------------------------------------------------------
898: dr=r0/(fa%r1-1)
899: dt1=twopi/fa%p2
900: dt3=twopi/(3*(fa%p1-fa%p2-1))
901: CALL barry_get_global_corners(fa,0,x,y,m,n)
902: !-----------------------------------------------------------------------
903: ! fill array.
904: !-----------------------------------------------------------------------
905: CALL VecGetArray(g,ga,ig,ierr)
906: k=ig+2*fa%offg(0)
907: IF(diagnose)THEN
908: WRITE(*,'(a)')"region 1"
909: WRITE(*,10)
910: ENDIF
911: DO j=y,y+n-1
912: r=2*r0+j*dr
913: DO i=x,x+m-1
914: IF(i < (fa%p1-fa%p2)/2)THEN
915: theta=i*dt3
916: ga(k+1)=r*COS(theta)
917: ga(k+2)=r*SIN(theta)-4*r0
918: ELSEIF(i > fa%p2 + (fa%p1-fa%p2)/2)then
919: theta=pi+i*dt3
920: ga(k+1)=r*COS(PETSC_PI+i*dt3)
921: ga(k+2)=r*SIN(PETSC_PI+i*dt3)-4*r0
922: ELSE
923: theta=(i-(fa%p1-fa%p2)/2)*dt1-pi/2
924: ga(k+1)=r*COS(theta)
925: ga(k+2)=r*SIN(theta)
926: ENDIF
927: IF(diagnose)WRITE(*,20)j,i,r,theta,ga(k+1),ga(k+2)
928: k=k+2
929: ENDDO
930: IF(diagnose)WRITE(*,10)
931: ENDDO
932: CALL VecRestoreArray(g,ga,ig,ierr)
933: !-----------------------------------------------------------------------
934: ! terminate.
935: !-----------------------------------------------------------------------
936: RETURN
937: END SUBROUTINE barry_map_region1
938: END MODULE barry_mod
939: !-----------------------------------------------------------------------
940: ! subprogram 12. barry_main.
941: ! main program.
942: !-----------------------------------------------------------------------
943: !-----------------------------------------------------------------------
944: ! declarations.
945: !-----------------------------------------------------------------------
946: PROGRAM barry_main
947: USE barry_mod
948: IMPLICIT NONE
950: TYPE(fa_type) :: fa
951: Vec :: g,l
952: !-----------------------------------------------------------------------
953: ! initialize and create structure, and get vectors.
954: !-----------------------------------------------------------------------
955: CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
956: CALL barry_create_fa(fa)
957: CALL barry_get_global_vector(fa,g)
958: CALL barry_get_local_vector(fa,l)
959: !-----------------------------------------------------------------------
960: ! map regions.
961: !-----------------------------------------------------------------------
962: CALL barry_map_region1(fa,g)
963: CALL barry_map_region2(fa,g)
964: CALL barry_map_region3(fa,g)
965: !-----------------------------------------------------------------------
966: ! graphic output.
967: !-----------------------------------------------------------------------
968: CALL barry_global_to_local(fa,g,l)
969: CALL barry_draw_fa(fa,g)
970: CALL barry_draw_fa(fa,l)
971: !-----------------------------------------------------------------------
972: ! clean up and finalize.
973: !-----------------------------------------------------------------------
974: CALL VecDestroy(g,ierr)
975: CALL VecDestroy(l,ierr)
976: CALL barry_destroy_fa(fa)
977: CALL PetscFinalize(ierr)
978: !-----------------------------------------------------------------------
979: ! terminate.
980: !-----------------------------------------------------------------------
981: END PROGRAM barry_main