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