Actual source code: ex6.c
2: static char help[] = "Demonstrates using 3 DA's to manage a slightly non-trivial grid";
4: #include petscda.h
5: #include petscsys.h
7: struct _p_FA {
8: MPI_Comm comm[3];
9: PetscInt xl[3],yl[3],ml[3],nl[3]; /* corners and sizes of local vector in DA */
10: PetscInt xg[3],yg[3],mg[3],ng[3]; /* corners and sizes of global vector in DA */
11: PetscInt offl[3],offg[3]; /* offset in local and global vector of region 1, 2 and 3 portions */
12: Vec g,l;
13: VecScatter vscat;
14: PetscInt p1,p2,r1,r2,r1g,r2g,sw;
15: };
16: typedef struct _p_FA *FA;
18: typedef struct {
19: PetscScalar X;
20: PetscScalar Y;
21: } Field;
23: PetscErrorCode FAGetLocalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
24: {
26: if (fa->comm[j]) {
27: *x = fa->xl[j];
28: *y = fa->yl[j];
29: *m = fa->ml[j];
30: *n = fa->nl[j];
31: } else {
32: *x = *y = *m = *n = 0;
33: }
34: return(0);
35: }
37: PetscErrorCode FAGetGlobalCorners(FA fa,PetscInt j,PetscInt *x,PetscInt *y,PetscInt *m,PetscInt *n)
38: {
40: if (fa->comm[j]) {
41: *x = fa->xg[j];
42: *y = fa->yg[j];
43: *m = fa->mg[j];
44: *n = fa->ng[j];
45: } else {
46: *x = *y = *m = *n = 0;
47: }
48: return(0);
49: }
51: PetscErrorCode FAGetLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
52: {
54: PetscScalar *va;
55: PetscInt i;
56: Field **a;
59: if (fa->comm[j]) {
60: VecGetArray(v,&va);
61: PetscMalloc(fa->nl[j]*sizeof(Field*),&a);
62: for (i=0; i<fa->nl[j]; i++) (a)[i] = (Field*) (va + 2*fa->offl[j] + i*2*fa->ml[j] - 2*fa->xl[j]);
63: *f = a - fa->yl[j];
64: VecRestoreArray(v,&va);
65: } else {
66: *f = 0;
67: }
68: return(0);
69: }
71: PetscErrorCode FARestoreLocalArray(FA fa,Vec v,PetscInt j,Field ***f)
72: {
74: void *dummy;
77: if (fa->comm[j]) {
78: dummy = *f + fa->yl[j];
79: PetscFree(dummy);
80: }
81: return(0);
82: }
84: PetscErrorCode FAGetGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
85: {
87: PetscScalar *va;
88: PetscInt i;
89: Field **a;
92: if (fa->comm[j]) {
93: VecGetArray(v,&va);
94: PetscMalloc(fa->ng[j]*sizeof(Field*),&a);
95: for (i=0; i<fa->ng[j]; i++) (a)[i] = (Field*) (va + 2*fa->offg[j] + i*2*fa->mg[j] - 2*fa->xg[j]);
96: *f = a - fa->yg[j];
97: VecRestoreArray(v,&va);
98: } else {
99: *f = 0;
100: }
101: return(0);
102: }
104: PetscErrorCode FARestoreGlobalArray(FA fa,Vec v,PetscInt j,Field ***f)
105: {
107: void *dummy;
110: if (fa->comm[j]) {
111: dummy = *f + fa->yg[j];
112: PetscFree(dummy);
113: }
114: return(0);
115: }
117: PetscErrorCode FAGetGlobalVector(FA fa,Vec *v)
118: {
121: VecDuplicate(fa->g,v);
122: return(0);
123: }
125: PetscErrorCode FAGetLocalVector(FA fa,Vec *v)
126: {
129: VecDuplicate(fa->l,v);
130: return(0);
131: }
133: PetscErrorCode FAGlobalToLocal(FA fa,Vec g,Vec l)
134: {
137: VecScatterBegin(g,l,INSERT_VALUES,SCATTER_FORWARD,fa->vscat);
138: VecScatterEnd(g,l,INSERT_VALUES,SCATTER_FORWARD,fa->vscat);
139: return(0);
140: }
142: PetscErrorCode FADestroy(FA fa)
143: {
146: VecDestroy(fa->g);
147: VecDestroy(fa->l);
148: VecScatterDestroy(fa->vscat);
149: PetscFree(fa);
150: return(0);
151: }
153: PetscErrorCode FACreate(FA *infa)
154: {
155: FA fa;
156: PetscMPIInt rank;
157: PetscInt tonglobal,globalrstart,x,nx,y,ny,*tonatural,i,j,*to,*from,offt[3];
158: PetscInt *fromnatural,fromnglobal,nscat,nlocal,cntl1,cntl2,cntl3,*indices;
161: /* Each DA manages the local vector for the portion of region 1, 2, and 3 for that processor
162: Each DA can belong on any subset (overlapping between DA's or not) of processors
163: For processes that a particular DA does not exist on, the corresponding comm should be set to zero
164: */
165: DA da1 = 0,da2 = 0,da3 = 0;
166: /*
167: v1, v2, v3 represent the local vector for a single DA
168: */
169: Vec vl1 = 0,vl2 = 0,vl3 = 0, vg1 = 0, vg2 = 0,vg3 = 0;
171: /*
172: globalvec and friends represent the global vectors that are used for the PETSc solvers
173: localvec represents the concatenation of the (up to) 3 local vectors; vl1, vl2, vl3
175: tovec and friends represent intermediate vectors that are ONLY used for setting up the
176: final communication patterns. Once this setup routine is complete they are destroyed.
177: The tovec is like the globalvec EXCEPT it has redundant locations for the ghost points
178: between regions 2+3 and 1.
179: */
180: AO toao,globalao;
181: IS tois,globalis,is;
182: Vec tovec,globalvec,localvec;
183: VecScatter vscat;
184: PetscScalar *globalarray,*localarray,*toarray;
186: PetscNew(struct _p_FA,&fa);
187: /*
188: fa->sw is the stencil width
190: fa->p1 is the width of region 1, fa->p2 the width of region 2 (must be the same)
191: fa->r1 height of region 1
192: fa->r2 height of region 2
193:
194: fa->r2 is also the height of region 3-4
195: (fa->p1 - fa->p2)/2 is the width of both region 3 and region 4
196: */
197: fa->p1 = 24;
198: fa->p2 = 15;
199: fa->r1 = 6;
200: fa->r2 = 6;
201: fa->sw = 1;
202: fa->r1g = fa->r1 + fa->sw;
203: fa->r2g = fa->r2 + fa->sw;
205: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
207: fa->comm[0] = PETSC_COMM_WORLD;
208: fa->comm[1] = PETSC_COMM_WORLD;
209: fa->comm[2] = PETSC_COMM_WORLD;
210: /* Test case with different communicators */
211: /* Normally one would use MPI_Comm routines to build MPI communicators on which you wish to partition the DAs*/
212: /*
213: if (rank == 0) {
214: fa->comm[0] = PETSC_COMM_SELF;
215: fa->comm[1] = 0;
216: fa->comm[2] = 0;
217: } else if (rank == 1) {
218: fa->comm[0] = 0;
219: fa->comm[1] = PETSC_COMM_SELF;
220: fa->comm[2] = 0;
221: } else {
222: fa->comm[0] = 0;
223: fa->comm[1] = 0;
224: fa->comm[2] = PETSC_COMM_SELF;
225: } */
227: if (fa->p2 > fa->p1 - 3) SETERRQ(1,"Width of region fa->p2 must be at least 3 less then width of region 1");
228: if (!((fa->p2 - fa->p1) % 2)) SETERRQ(1,"width of region 3 must NOT be divisible by 2");
230: if (fa->comm[1]) {
231: DACreate2d(fa->comm[1],DA_XPERIODIC,DA_STENCIL_BOX,fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da2);
232: DAGetLocalVector(da2,&vl2);
233: DAGetGlobalVector(da2,&vg2);
234: }
235: if (fa->comm[2]) {
236: DACreate2d(fa->comm[2],DA_NONPERIODIC,DA_STENCIL_BOX,fa->p1-fa->p2,fa->r2g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da3);
237: DAGetLocalVector(da3,&vl3);
238: DAGetGlobalVector(da3,&vg3);
239: }
240: if (fa->comm[0]) {
241: DACreate2d(fa->comm[0],DA_NONPERIODIC,DA_STENCIL_BOX,fa->p1,fa->r1g,PETSC_DECIDE,PETSC_DECIDE,1,fa->sw,PETSC_NULL,PETSC_NULL,&da1);
242: DAGetLocalVector(da1,&vl1);
243: DAGetGlobalVector(da1,&vg1);
244: }
246: /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership
247: for global vector with redundancy */
248: tonglobal = 0;
249: if (fa->comm[1]) {
250: DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
251: tonglobal += nx*ny;
252: }
253: if (fa->comm[2]) {
254: DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
255: tonglobal += nx*ny;
256: }
257: if (fa->comm[0]) {
258: DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
259: tonglobal += nx*ny;
260: }
261: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,tonglobal);
262: PetscSynchronizedFlush(PETSC_COMM_WORLD);
263:
264: /* Get tonatural number for each node */
265: PetscMalloc((tonglobal+1)*sizeof(PetscInt),&tonatural);
266: tonglobal = 0;
267: if (fa->comm[1]) {
268: DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
269: for (j=0; j<ny; j++) {
270: for (i=0; i<nx; i++) {
271: tonatural[tonglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
272: }
273: }
274: }
275: if (fa->comm[2]) {
276: DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
277: for (j=0; j<ny; j++) {
278: for (i=0; i<nx; i++) {
279: if (x + i < (fa->p1 - fa->p2)/2) tonatural[tonglobal++] = x + i + fa->p1*(y + j);
280: else tonatural[tonglobal++] = fa->p2 + x + i + fa->p1*(y + j);
281: }
282: }
283: }
284: if (fa->comm[0]) {
285: DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
286: for (j=0; j<ny; j++) {
287: for (i=0; i<nx; i++) {
288: tonatural[tonglobal++] = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
289: }
290: }
291: }
292: /* PetscIntView(tonglobal,tonatural,PETSC_VIEWER_STDOUT_WORLD); */
293: AOCreateBasic(PETSC_COMM_WORLD,tonglobal,tonatural,0,&toao);
294: PetscFree(tonatural);
296: /* count the number of unknowns owned on each processor and determine the starting point of each processors ownership
297: for global vector without redundancy */
298: fromnglobal = 0;
299: fa->offg[1] = 0;
300: offt[1] = 0;
301: if (fa->comm[1]) {
302: DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
303: offt[2] = nx*ny;
304: if (y+ny == fa->r2g) {ny--;} /* includes the ghost points on the upper side */
305: fromnglobal += nx*ny;
306: fa->offg[2] = fromnglobal;
307: } else {
308: offt[2] = 0;
309: fa->offg[2] = 0;
310: }
311: if (fa->comm[2]) {
312: DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
313: offt[0] = offt[2] + nx*ny;
314: if (y+ny == fa->r2g) {ny--;} /* includes the ghost points on the upper side */
315: fromnglobal += nx*ny;
316: fa->offg[0] = fromnglobal;
317: } else {
318: offt[0] = offt[2];
319: fa->offg[0] = fromnglobal;
320: }
321: if (fa->comm[0]) {
322: DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
323: if (y == 0) {ny--;} /* includes the ghost points on the lower side */
324: fromnglobal += nx*ny;
325: }
326: MPI_Scan(&fromnglobal,&globalrstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
327: globalrstart -= fromnglobal;
328: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Number of unknowns owned %d\n",rank,fromnglobal);
329: PetscSynchronizedFlush(PETSC_COMM_WORLD);
331: /* Get fromnatural number for each node */
332: PetscMalloc((fromnglobal+1)*sizeof(PetscInt),&fromnatural);
333: fromnglobal = 0;
334: if (fa->comm[1]) {
335: DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
336: if (y+ny == fa->r2g) {ny--;} /* includes the ghost points on the upper side */
337: fa->xg[1] = x; fa->yg[1] = y; fa->mg[1] = nx; fa->ng[1] = ny;
338: DAGetGhostCorners(da2,&fa->xl[1],&fa->yl[1],0,&fa->ml[1],&fa->nl[1],0);
339: for (j=0; j<ny; j++) {
340: for (i=0; i<nx; i++) {
341: fromnatural[fromnglobal++] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);
342: }
343: }
344: }
345: if (fa->comm[2]) {
346: DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
347: if (y+ny == fa->r2g) {ny--;} /* includes the ghost points on the upper side */
348: fa->xg[2] = x; fa->yg[2] = y; fa->mg[2] = nx; fa->ng[2] = ny;
349: DAGetGhostCorners(da3,&fa->xl[2],&fa->yl[2],0,&fa->ml[2],&fa->nl[2],0);
350: for (j=0; j<ny; j++) {
351: for (i=0; i<nx; i++) {
352: if (x + i < (fa->p1 - fa->p2)/2) fromnatural[fromnglobal++] = x + i + fa->p1*(y + j);
353: else fromnatural[fromnglobal++] = fa->p2 + x + i + fa->p1*(y + j);
354: }
355: }
356: }
357: if (fa->comm[0]) {
358: DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
359: if (y == 0) {ny--;} /* includes the ghost points on the lower side */
360: else y--;
361: fa->xg[0] = x; fa->yg[0] = y; fa->mg[0] = nx; fa->ng[0] = ny;
362: DAGetGhostCorners(da1,&fa->xl[0],&fa->yl[0],0,&fa->ml[0],&fa->nl[0],0);
363: for (j=0; j<ny; j++) {
364: for (i=0; i<nx; i++) {
365: fromnatural[fromnglobal++] = fa->p1*fa->r2 + x + i + fa->p1*(y + j);
366: }
367: }
368: }
369: /*PetscIntView(fromnglobal,fromnatural,PETSC_VIEWER_STDOUT_WORLD);*/
370: AOCreateBasic(PETSC_COMM_WORLD,fromnglobal,fromnatural,0,&globalao);
371: PetscFree(fromnatural);
373: /* ---------------------------------------------------*/
374: /* Create the scatter that updates 1 from 2 and 3 and 3 and 2 from 1 */
375: /* currently handles stencil width of 1 ONLY */
376: PetscMalloc(tonglobal*sizeof(PetscInt),&to);
377: PetscMalloc(tonglobal*sizeof(PetscInt),&from);
378: nscat = 0;
379: if (fa->comm[1]) {
380: DAGetCorners(da2,&x,&y,0,&nx,&ny,0);
381: for (j=0; j<ny; j++) {
382: for (i=0; i<nx; i++) {
383: to[nscat] = from[nscat] = (fa->p1 - fa->p2)/2 + x + i + fa->p1*(y + j);nscat++;
384: }
385: }
386: }
387: if (fa->comm[2]) {
388: DAGetCorners(da3,&x,&y,0,&nx,&ny,0);
389: for (j=0; j<ny; j++) {
390: for (i=0; i<nx; i++) {
391: if (x + i < (fa->p1 - fa->p2)/2) {
392: to[nscat] = from[nscat] = x + i + fa->p1*(y + j);nscat++;
393: } else {
394: to[nscat] = from[nscat] = fa->p2 + x + i + fa->p1*(y + j);nscat++;
395: }
396: }
397: }
398: }
399: if (fa->comm[0]) {
400: DAGetCorners(da1,&x,&y,0,&nx,&ny,0);
401: for (j=0; j<ny; j++) {
402: for (i=0; i<nx; i++) {
403: to[nscat] = fa->p1*fa->r2g + x + i + fa->p1*(y + j);
404: from[nscat++] = fa->p1*(fa->r2 - 1) + x + i + fa->p1*(y + j);
405: }
406: }
407: }
408: AOApplicationToPetsc(toao,nscat,to);
409: AOApplicationToPetsc(globalao,nscat,from);
410: ISCreateGeneral(PETSC_COMM_WORLD,nscat,to,&tois);
411: ISCreateGeneral(PETSC_COMM_WORLD,nscat,from,&globalis);
412: PetscFree(to);
413: PetscFree(from);
414: VecCreateMPI(PETSC_COMM_WORLD,tonglobal,PETSC_DETERMINE,&tovec);
415: VecCreateMPI(PETSC_COMM_WORLD,fromnglobal,PETSC_DETERMINE,&globalvec);
416: VecScatterCreate(globalvec,globalis,tovec,tois,&vscat);
417: ISDestroy(tois);
418: ISDestroy(globalis);
419: AODestroy(globalao);
420: AODestroy(toao);
422: /* fill up global vector without redundant values with PETSc global numbering */
423: VecGetArray(globalvec,&globalarray);
424: for (i=0; i<fromnglobal; i++) {
425: globalarray[i] = globalrstart + i;
426: }
427: VecRestoreArray(globalvec,&globalarray);
428:
429: /* scatter PETSc global indices to redundant valueed array */
430: VecScatterBegin(globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD,vscat);
431: VecScatterEnd(globalvec,tovec,INSERT_VALUES,SCATTER_FORWARD,vscat);
432:
433: /* Create local vector that is the concatenation of the local vectors */
434: nlocal = 0;
435: cntl1 = cntl2 = cntl3 = 0;
436: if (fa->comm[1]) {
437: VecGetSize(vl2,&cntl2);
438: nlocal += cntl2;
439: }
440: if (fa->comm[2]) {
441: VecGetSize(vl3,&cntl3);
442: nlocal += cntl3;
443: }
444: if (fa->comm[0]) {
445: VecGetSize(vl1,&cntl1);
446: nlocal += cntl1;
447: }
448: fa->offl[0] = cntl2 + cntl3;
449: fa->offl[1] = 0;
450: fa->offl[2] = cntl2;
451: VecCreateSeq(PETSC_COMM_SELF,nlocal,&localvec);
452:
453: /* cheat so that vl1, vl2, vl3 shared array memory with localvec */
454: VecGetArray(localvec,&localarray);
455: VecGetArray(tovec,&toarray);
456: if (fa->comm[1]) {
457: VecPlaceArray(vl2,localarray+fa->offl[1]);
458: VecPlaceArray(vg2,toarray+offt[1]);
459: DAGlobalToLocalBegin(da2,vg2,INSERT_VALUES,vl2);
460: DAGlobalToLocalEnd(da2,vg2,INSERT_VALUES,vl2);
461: DARestoreGlobalVector(da2,&vg2);
462: }
463: if (fa->comm[2]) {
464: VecPlaceArray(vl3,localarray+fa->offl[2]);
465: VecPlaceArray(vg3,toarray+offt[2]);
466: DAGlobalToLocalBegin(da3,vg3,INSERT_VALUES,vl3);
467: DAGlobalToLocalEnd(da3,vg3,INSERT_VALUES,vl3);
468: DARestoreGlobalVector(da3,&vg3);
469: }
470: if (fa->comm[0]) {
471: VecPlaceArray(vl1,localarray+fa->offl[0]);
472: VecPlaceArray(vg1,toarray+offt[0]);
473: DAGlobalToLocalBegin(da1,vg1,INSERT_VALUES,vl1);
474: DAGlobalToLocalEnd(da1,vg1,INSERT_VALUES,vl1);
475: DARestoreGlobalVector(da1,&vg1);
476: }
477: VecRestoreArray(localvec,&localarray);
478: VecRestoreArray(tovec,&toarray);
480: /* no longer need the redundant vector and VecScatter to it */
481: VecScatterDestroy(vscat);
482: VecDestroy(tovec);
484: /* Create final scatter that goes directly from globalvec to localvec */
485: /* this is the one to be used in the application code */
486: PetscMalloc(nlocal*sizeof(PetscInt),&indices);
487: VecGetArray(localvec,&localarray);
488: for (i=0; i<nlocal; i++) {
489: indices[i] = (PetscInt) (2*localarray[i]);
490: }
491: VecRestoreArray(localvec,&localarray);
492: ISCreateBlock(PETSC_COMM_WORLD,2,nlocal,indices,&is);
493: PetscFree(indices);
495: VecCreateSeq(PETSC_COMM_SELF,2*nlocal,&fa->l);
496: VecCreateMPI(PETSC_COMM_WORLD,2*fromnglobal,PETSC_DETERMINE,&fa->g);
498: VecScatterCreate(fa->g,is,fa->l,PETSC_NULL,&fa->vscat);
499: ISDestroy(is);
501: VecDestroy(globalvec);
502: VecDestroy(localvec);
503: if (fa->comm[0]) {
504: DARestoreLocalVector(da1,&vl1);
505: DADestroy(da1);
506: }
507: if (fa->comm[1]) {
508: DARestoreLocalVector(da2,&vl2);
509: DADestroy(da2);
510: }
511: if (fa->comm[2]) {
512: DARestoreLocalVector(da3,&vl3);
513: DADestroy(da3);
514: }
515: *infa = fa;
516: return(0);
517: }
519: /* Crude graphics to test that the ghost points are properly updated */
520: #include petscdraw.h
522: typedef struct {
523: PetscInt m[3],n[3];
524: PetscScalar *xy[3];
525: } ZoomCtx;
527: PetscErrorCode DrawPatch(PetscDraw draw,void *ctx)
528: {
529: ZoomCtx *zctx = (ZoomCtx*)ctx;
531: PetscInt m,n,i,j,id,k;
532: PetscReal x1,x2,x3,x4,y_1,y2,y3,y4;
533: PetscScalar *xy;
536: for (k=0; k<3; k++) {
537: m = zctx->m[k];
538: n = zctx->n[k];
539: xy = zctx->xy[k];
541: for (j=0; j<n-1; j++) {
542: for (i=0; i<m-1; i++) {
543: id = i+j*m; x1 = xy[2*id];y_1 = xy[2*id+1];
544: id = i+j*m+1; x2 = xy[2*id];y2 = xy[2*id+1];
545: id = i+j*m+1+m;x3 = xy[2*id];y3 = xy[2*id+1];
546: id = i+j*m+m; x4 = xy[2*id];y4 = xy[2*id+1];
547: PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);
548: PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);
549: PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);
550: PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);
551: }
552: }
553: }
554: PetscDrawFlush(draw);
555: return(0);
556: }
558: PetscErrorCode DrawFA(FA fa,Vec v)
559: {
561: PetscScalar *va;
562: ZoomCtx zctx;
563: PetscDraw draw;
564: PetscReal xmint = 10000.0,xmaxt = -10000.0,ymint = 100000.0,ymaxt = -10000.0;
565: PetscReal xmin,xmax,ymin,ymax;
566: PetscInt i,vn,ln,j;
569: VecGetArray(v,&va);
570: VecGetSize(v,&vn);
571: VecGetSize(fa->l,&ln);
572: for (j=0; j<3; j++) {
573: if (vn == ln) {
574: zctx.xy[j] = va + 2*fa->offl[j];
575: zctx.m[j] = fa->ml[j];
576: zctx.n[j] = fa->nl[j];
577: } else {
578: zctx.xy[j] = va + 2*fa->offg[j];
579: zctx.m[j] = fa->mg[j];
580: zctx.n[j] = fa->ng[j];
581: }
582: for (i=0; i<zctx.m[j]*zctx.n[j]; i++) {
583: if (zctx.xy[j][2*i] > xmax) xmax = zctx.xy[j][2*i];
584: if (zctx.xy[j][2*i] < xmin) xmin = zctx.xy[j][2*i];
585: if (zctx.xy[j][2*i+1] > ymax) ymax = zctx.xy[j][2*i+1];
586: if (zctx.xy[j][2*i+1] < ymin) ymin = zctx.xy[j][2*i+1];
587: }
588: }
589: MPI_Allreduce(&xmin,&xmint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
590: MPI_Allreduce(&xmax,&xmaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
591: MPI_Allreduce(&ymin,&ymint,1,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD);
592: MPI_Allreduce(&ymax,&ymaxt,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
593: xmin = xmint - .2*(xmaxt - xmint);
594: xmax = xmaxt + .2*(xmaxt - xmint);
595: ymin = ymint - .2*(ymaxt - ymint);
596: ymax = ymaxt + .2*(ymaxt - ymint);
597: PetscDrawOpenX(PETSC_COMM_WORLD,0,"meshes",PETSC_DECIDE,PETSC_DECIDE,700,700,&draw);
598: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
599: PetscDrawZoom(draw,DrawPatch,&zctx);
600: VecRestoreArray(v,&va);
601: PetscDrawDestroy(draw);
602: return(0);
603: }
605: /* crude mappings from rectangular arrays to the true geometry. These are ONLY for testing!
606: they will not be used the actual code */
607: PetscErrorCode FAMapRegion3(FA fa,Vec g)
608: {
610: PetscReal R = 1.0,Rscale,Ascale;
611: PetscInt i,k,x,y,m,n;
612: Field **ga;
615: Rscale = R/(fa->r2-1);
616: Ascale = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));
618: FAGetGlobalCorners(fa,2,&x,&y,&m,&n);
619: FAGetGlobalArray(fa,g,2,&ga);
620: for (k=y; k<y+n; k++) {
621: for (i=x; i<x+m; i++) {
622: ga[k][i].X = (R + k*Rscale)*PetscCosScalar(1.*PETSC_PI/6. + i*Ascale);
623: ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(1.*PETSC_PI/6. + i*Ascale) - 4.*R;
624: }
625: }
626: FARestoreGlobalArray(fa,g,2,&ga);
627: return(0);
628: }
630: PetscErrorCode FAMapRegion2(FA fa,Vec g)
631: {
633: PetscReal R = 1.0,Rscale,Ascale;
634: PetscInt i,k,x,y,m,n;
635: Field **ga;
638: Rscale = R/(fa->r2-1);
639: Ascale = 2.0*PETSC_PI/fa->p2;
641: FAGetGlobalCorners(fa,1,&x,&y,&m,&n);
642: FAGetGlobalArray(fa,g,1,&ga);
643: for (k=y; k<y+n; k++) {
644: for (i=x; i<x+m; i++) {
645: ga[k][i].X = (R + k*Rscale)*PetscCosScalar(i*Ascale - PETSC_PI/2.0);
646: ga[k][i].Y = (R + k*Rscale)*PetscSinScalar(i*Ascale - PETSC_PI/2.0);
647: }
648: }
649: FARestoreGlobalArray(fa,g,1,&ga);
650: return(0);
651: }
653: PetscErrorCode FAMapRegion1(FA fa,Vec g)
654: {
656: PetscReal R = 1.0,Rscale,Ascale1,Ascale3;
657: PetscInt i,k,x,y,m,n;
658: Field **ga;
661: Rscale = R/(fa->r1-1);
662: Ascale1 = 2.0*PETSC_PI/fa->p2;
663: Ascale3 = 2.0*PETSC_PI/(3.0*(fa->p1 - fa->p2 - 1));
665: FAGetGlobalCorners(fa,0,&x,&y,&m,&n);
666: FAGetGlobalArray(fa,g,0,&ga);
668: /* This mapping is WRONG! Not sure how to do it so I've done a poor job of
669: it. You can see that the grid connections are correct. */
670: for (k=y; k<y+n; k++) {
671: for (i=x; i<x+m; i++) {
672: if (i < (fa->p1-fa->p2)/2) {
673: ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(i*Ascale3);
674: ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(i*Ascale3) - 4.*R;
675: } else if (i > fa->p2 + (fa->p1 - fa->p2)/2) {
676: ga[k][i].X = (2.0*R + k*Rscale)*PetscCosScalar(PETSC_PI+i*Ascale3);
677: ga[k][i].Y = (2.0*R + k*Rscale)*PetscSinScalar(PETSC_PI+i*Ascale3) - 4.*R;
678: } else {
679: ga[k][i].X = (2.*R + k*Rscale)*PetscCosScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
680: ga[k][i].Y = (2.*R + k*Rscale)*PetscSinScalar((i-(fa->p1-fa->p2)/2)*Ascale1 - PETSC_PI/2.0);
681: }
682: }
683: }
684: FARestoreGlobalArray(fa,g,0,&ga);
685: return(0);
686: }
688: /* Simple test to check that the ghost points are properly updated */
689: PetscErrorCode FATest(FA fa)
690: {
692: Vec l,g;
693: Field **la;
694: PetscInt x,y,m,n,j,i,k,p;
695: PetscMPIInt rank;
698: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
700: FAGetGlobalVector(fa,&g);
701: FAGetLocalVector(fa,&l);
703: /* fill up global vector of one region at a time with ITS logical coordinates, then update LOCAL
704: vector; print local vectors to confirm they are correctly filled */
705: for (j=0; j<3; j++) {
706: VecSet(g,0.0);
707: FAGetGlobalCorners(fa,j,&x,&y,&m,&n);
708: PetscPrintf(PETSC_COMM_WORLD,"\nFilling global region %d, showing local results \n",j+1);
709: FAGetGlobalArray(fa,g,j,&la);
710: for (k=y; k<y+n; k++) {
711: for (i=x; i<x+m; i++) {
712: la[k][i].X = i;
713: la[k][i].Y = k;
714: }
715: }
716: FARestoreGlobalArray(fa,g,j,&la);
718: FAGlobalToLocal(fa,g,l);
719: DrawFA(fa,g);
720: DrawFA(fa,l);
722: for (p=0; p<3; p++) {
723: FAGetLocalCorners(fa,p,&x,&y,&m,&n);
724: FAGetLocalArray(fa,l,p,&la);
725: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n[%d] Local array for region %d \n",rank,p+1);
726: for (k=y+n-1; k>=y; k--) { /* print in reverse order to match diagram in paper */
727: for (i=x; i<x+m; i++) {
728: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"(%G,%G) ",la[k][i].X,la[k][i].Y);
729: }
730: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"\n");
731: }
732: FARestoreLocalArray(fa,l,p,&la);
733: PetscSynchronizedFlush(PETSC_COMM_WORLD);
734: }
735: }
736: VecDestroy(g);
737: VecDestroy(l);
738: return(0);
739: }
743: int main(int argc,char **argv)
744: {
745: FA fa;
747: Vec g,l;
749: PetscInitialize(&argc,&argv,0,help);
751: FACreate(&fa);
752: /* FATest(fa);*/
754: FAGetGlobalVector(fa,&g);
755: FAGetLocalVector(fa,&l);
757: FAMapRegion1(fa,g);
758: FAMapRegion2(fa,g);
759: FAMapRegion3(fa,g);
761: FAGlobalToLocal(fa,g,l);
762: DrawFA(fa,g);
763: DrawFA(fa,l);
765: VecDestroy(g);
766: VecDestroy(l);
767: FADestroy(fa);
769: PetscFinalize();
770: return 0;
771: }