Actual source code: da2.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/
6: PetscErrorCode DMView_DA_2d(DM da,PetscViewer viewer)
7: {
9: PetscMPIInt rank;
10: PetscBool iascii,isdraw,isbinary;
11: DM_DA *dd = (DM_DA*)da->data;
12: #if defined(PETSC_HAVE_MATLAB_ENGINE)
13: PetscBool ismatlab;
14: #endif
17: MPI_Comm_rank(((PetscObject)da)->comm,&rank);
19: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
20: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
21: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
22: #if defined(PETSC_HAVE_MATLAB_ENGINE)
23: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
24: #endif
25: if (iascii) {
26: PetscViewerFormat format;
28: PetscViewerGetFormat(viewer, &format);
29: if (format != PETSC_VIEWER_ASCII_VTK && format != PETSC_VIEWER_ASCII_VTK_CELL) {
30: DMDALocalInfo info;
31: DMDAGetLocalInfo(da,&info);
32: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
33: PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,dd->M,dd->N,dd->m,dd->n,dd->w,dd->s);
34: PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",info.xs,info.xs+info.xm,info.ys,info.ys+info.ym);
35: PetscViewerFlush(viewer);
36: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
37: } else {
38: DMView_DA_VTK(da,viewer);
39: }
40: } else if (isdraw) {
41: PetscDraw draw;
42: double ymin = -1*dd->s-1,ymax = dd->N+dd->s;
43: double xmin = -1*dd->s-1,xmax = dd->M+dd->s;
44: double x,y;
45: PetscInt base,*idx;
46: char node[10];
47: PetscBool isnull;
48:
49: PetscViewerDrawGetDraw(viewer,0,&draw);
50: PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
51: if (!dd->coordinates) {
52: PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
53: }
54: PetscDrawSynchronizedClear(draw);
56: /* first processor draw all node lines */
57: if (!rank) {
58: ymin = 0.0; ymax = dd->N - 1;
59: for (xmin=0; xmin<dd->M; xmin++) {
60: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
61: }
62: xmin = 0.0; xmax = dd->M - 1;
63: for (ymin=0; ymin<dd->N; ymin++) {
64: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
65: }
66: }
67: PetscDrawSynchronizedFlush(draw);
68: PetscDrawPause(draw);
70: /* draw my box */
71: ymin = dd->ys; ymax = dd->ye - 1; xmin = dd->xs/dd->w;
72: xmax =(dd->xe-1)/dd->w;
73: PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
74: PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
75: PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
76: PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);
78: /* put in numbers */
79: base = (dd->base)/dd->w;
80: for (y=ymin; y<=ymax; y++) {
81: for (x=xmin; x<=xmax; x++) {
82: sprintf(node,"%d",(int)base++);
83: PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
84: }
85: }
87: PetscDrawSynchronizedFlush(draw);
88: PetscDrawPause(draw);
89: /* overlay ghost numbers, useful for error checking */
90: /* put in numbers */
92: base = 0; idx = dd->idx;
93: ymin = dd->Ys; ymax = dd->Ye; xmin = dd->Xs; xmax = dd->Xe;
94: for (y=ymin; y<ymax; y++) {
95: for (x=xmin; x<xmax; x++) {
96: if ((base % dd->w) == 0) {
97: sprintf(node,"%d",(int)(idx[base]/dd->w));
98: PetscDrawString(draw,x/dd->w,y,PETSC_DRAW_BLUE,node);
99: }
100: base++;
101: }
102: }
103: PetscDrawSynchronizedFlush(draw);
104: PetscDrawPause(draw);
105: } else if (isbinary){
106: DMView_DA_Binary(da,viewer);
107: #if defined(PETSC_HAVE_MATLAB_ENGINE)
108: } else if (ismatlab) {
109: DMView_DA_Matlab(da,viewer);
110: #endif
111: } else SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for DMDA 1d",((PetscObject)viewer)->type_name);
112: return(0);
113: }
115: /*
116: M is number of grid points
117: m is number of processors
119: */
122: PetscErrorCode DMDASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
123: {
125: PetscInt m,n = 0,x = 0,y = 0;
126: PetscMPIInt size,csize,rank;
129: MPI_Comm_size(comm,&size);
130: MPI_Comm_rank(comm,&rank);
132: csize = 4*size;
133: do {
134: if (csize % 4) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y);
135: csize = csize/4;
136:
137: m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N)));
138: if (!m) m = 1;
139: while (m > 0) {
140: n = csize/m;
141: if (m*n == csize) break;
142: m--;
143: }
144: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
146: x = M/m + ((M % m) > ((csize-1) % m));
147: y = (N + (csize-1)/m)/n;
148: } while ((x < 4 || y < 4) && csize > 1);
149: if (size != csize) {
150: MPI_Group entire_group,sub_group;
151: PetscMPIInt i,*groupies;
153: MPI_Comm_group(comm,&entire_group);
154: PetscMalloc(csize*sizeof(PetscInt),&groupies);
155: for (i=0; i<csize; i++) {
156: groupies[i] = (rank/csize)*csize + i;
157: }
158: MPI_Group_incl(entire_group,csize,groupies,&sub_group);
159: PetscFree(groupies);
160: MPI_Comm_create(comm,sub_group,outcomm);
161: MPI_Group_free(&entire_group);
162: MPI_Group_free(&sub_group);
163: PetscInfo1(0,"DMDASplitComm2d:Creating redundant coarse problems of size %d\n",csize);
164: } else {
165: *outcomm = comm;
166: }
167: return(0);
168: }
172: static PetscErrorCode DMDAFunction(DM dm,Vec x,Vec F)
173: {
175: Vec localX;
176:
178: DMGetLocalVector(dm,&localX);
179: DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
180: DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
181: DMDAComputeFunction1(dm,localX,F,dm->ctx);
182: DMRestoreLocalVector(dm,&localX);
183: return(0);
184: }
188: /*@C
189: DMDASetLocalFunction - Caches in a DM a local function.
191: Logically Collective on DMDA
193: Input Parameter:
194: + da - initial distributed array
195: - lf - the local function
197: Level: intermediate
199: Notes:
200: If you use SNESSetDM() and did not use SNESSetFunction() then the context passed to your function is the context set with DMSetApplicationContext()
202: Developer Notes: It is possibly confusing which context is passed to the user function, it would be nice to unify them somehow.
204: .keywords: distributed array, refine
206: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunctioni()
207: @*/
208: PetscErrorCode DMDASetLocalFunction(DM da,DMDALocalFunction1 lf)
209: {
211: DM_DA *dd = (DM_DA*)da->data;
215: DMSetFunction(da,DMDAFunction);
216: dd->lf = lf;
217: return(0);
218: }
222: /*@C
223: DMDASetLocalFunctioni - Caches in a DM a local function that evaluates a single component
225: Logically Collective on DMDA
227: Input Parameter:
228: + da - initial distributed array
229: - lfi - the local function
231: Level: intermediate
233: .keywords: distributed array, refine
235: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
236: @*/
237: PetscErrorCode DMDASetLocalFunctioni(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
238: {
239: DM_DA *dd = (DM_DA*)da->data;
242: dd->lfi = lfi;
243: return(0);
244: }
248: /*@C
249: DMDASetLocalFunctionib - Caches in a DM a block local function that evaluates a single component
251: Logically Collective on DMDA
253: Input Parameter:
254: + da - initial distributed array
255: - lfi - the local function
257: Level: intermediate
259: .keywords: distributed array, refine
261: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
262: @*/
263: PetscErrorCode DMDASetLocalFunctionib(DM da,PetscErrorCode (*lfi)(DMDALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
264: {
265: DM_DA *dd = (DM_DA*)da->data;
268: dd->lfib = lfi;
269: return(0);
270: }
274: PetscErrorCode DMDASetLocalAdicFunction_Private(DM da,DMDALocalFunction1 ad_lf)
275: {
276: DM_DA *dd = (DM_DA*)da->data;
279: dd->adic_lf = ad_lf;
280: return(0);
281: }
283: /*MC
284: DMDASetLocalAdicFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR
286: Synopsis:
287: PetscErrorCode DMDASetLocalAdicFunctioni(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
288:
289: Logically Collective on DMDA
291: Input Parameter:
292: + da - initial distributed array
293: - ad_lfi - the local function as computed by ADIC/ADIFOR
295: Level: intermediate
297: .keywords: distributed array, refine
299: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
300: DMDASetLocalJacobian(), DMDASetLocalFunctioni()
301: M*/
305: PetscErrorCode DMDASetLocalAdicFunctioni_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
306: {
307: DM_DA *dd = (DM_DA*)da->data;
310: dd->adic_lfi = ad_lfi;
311: return(0);
312: }
314: /*MC
315: DMDASetLocalAdicMFFunctioni - Caches in a DM a local functioni computed by ADIC/ADIFOR
317: Synopsis:
318: PetscErrorCode DMDASetLocalAdicFunctioni(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
319:
320: Logically Collective on DMDA
322: Input Parameter:
323: + da - initial distributed array
324: - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR
326: Level: intermediate
328: .keywords: distributed array, refine
330: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
331: DMDASetLocalJacobian(), DMDASetLocalFunctioni()
332: M*/
336: PetscErrorCode DMDASetLocalAdicMFFunctioni_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
337: {
338: DM_DA *dd = (DM_DA*)da->data;
341: dd->adicmf_lfi = admf_lfi;
342: return(0);
343: }
345: /*MC
346: DMDASetLocalAdicFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR
348: Synopsis:
349: PetscErrorCode DMDASetLocalAdicFunctionib(DM da,PetscInt (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
350:
351: Logically Collective on DMDA
353: Input Parameter:
354: + da - initial distributed array
355: - ad_lfi - the local function as computed by ADIC/ADIFOR
357: Level: intermediate
359: .keywords: distributed array, refine
361: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
362: DMDASetLocalJacobian(), DMDASetLocalFunctionib()
363: M*/
367: PetscErrorCode DMDASetLocalAdicFunctionib_Private(DM da,PetscErrorCode (*ad_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
368: {
369: DM_DA *dd = (DM_DA*)da->data;
372: dd->adic_lfib = ad_lfi;
373: return(0);
374: }
376: /*MC
377: DMDASetLocalAdicMFFunctionib - Caches in a DM a block local functioni computed by ADIC/ADIFOR
379: Synopsis:
380: PetscErrorCode DMDASetLocalAdicFunctionib(DM da,int (ad_lf*)(DMDALocalInfo*,MatStencil*,void*,void*,void*)
382: Logically Collective on DMDA
384: Input Parameter:
385: + da - initial distributed array
386: - admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR
388: Level: intermediate
390: .keywords: distributed array, refine
392: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
393: DMDASetLocalJacobian(), DMDASetLocalFunctionib()
394: M*/
398: PetscErrorCode DMDASetLocalAdicMFFunctionib_Private(DM da,PetscErrorCode (*admf_lfi)(DMDALocalInfo*,MatStencil*,void*,void*,void*))
399: {
400: DM_DA *dd = (DM_DA*)da->data;
403: dd->adicmf_lfib = admf_lfi;
404: return(0);
405: }
407: /*MC
408: DMDASetLocalAdicMFFunction - Caches in a DM a local function computed by ADIC/ADIFOR
410: Synopsis:
411: PetscErrorCode DMDASetLocalAdicMFFunction(DM da,DMDALocalFunction1 ad_lf)
413: Logically Collective on DMDA
415: Input Parameter:
416: + da - initial distributed array
417: - ad_lf - the local function as computed by ADIC/ADIFOR
419: Level: intermediate
421: .keywords: distributed array, refine
423: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction(),
424: DMDASetLocalJacobian()
425: M*/
429: PetscErrorCode DMDASetLocalAdicMFFunction_Private(DM da,DMDALocalFunction1 ad_lf)
430: {
431: DM_DA *dd = (DM_DA*)da->data;
434: dd->adicmf_lf = ad_lf;
435: return(0);
436: }
440: PetscErrorCode DMDAJacobianLocal(DM dm,Vec x,Mat A,Mat B, MatStructure *str)
441: {
443: Vec localX;
444:
446: DMGetLocalVector(dm,&localX);
447: DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
448: DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
449: MatFDColoringApply(B,dm->fd,localX,str,dm);
450: DMRestoreLocalVector(dm,&localX);
451: /* Assemble true Jacobian; if it is different */
452: if (A != B) {
453: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
454: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
455: }
456: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
457: *str = SAME_NONZERO_PATTERN;
458: return(0);
459: }
464: static PetscErrorCode DMDAJacobian(DM dm,Vec x,Mat A,Mat B, MatStructure *str)
465: {
467: Vec localX;
468:
470: DMGetLocalVector(dm,&localX);
471: DMGlobalToLocalBegin(dm,x,INSERT_VALUES,localX);
472: DMGlobalToLocalEnd(dm,x,INSERT_VALUES,localX);
473: DMDAComputeJacobian1(dm,localX,B,dm->ctx);
474: DMRestoreLocalVector(dm,&localX);
475: /* Assemble true Jacobian; if it is different */
476: if (A != B) {
477: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
478: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
479: }
480: MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
481: *str = SAME_NONZERO_PATTERN;
482: return(0);
483: }
485: /*@C
486: DMDASetLocalJacobian - Caches in a DM a local Jacobian computation function
488: Logically Collective on DMDA
490:
491: Input Parameter:
492: + da - initial distributed array
493: - lj - the local Jacobian
495: Level: intermediate
497: .keywords: distributed array, refine
499: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalFunction()
500: @*/
503: PetscErrorCode DMDASetLocalJacobian(DM da,DMDALocalFunction1 lj)
504: {
506: DM_DA *dd = (DM_DA*)da->data;
510: DMSetJacobian(da,DMDAJacobian);
511: dd->lj = lj;
512: return(0);
513: }
517: /*@C
518: DMDAGetLocalFunction - Gets from a DM a local function and its ADIC/ADIFOR Jacobian
520: Note Collective
522: Input Parameter:
523: . da - initial distributed array
525: Output Parameter:
526: . lf - the local function
528: Level: intermediate
530: .keywords: distributed array, refine
532: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalJacobian(), DMDASetLocalFunction()
533: @*/
534: PetscErrorCode DMDAGetLocalFunction(DM da,DMDALocalFunction1 *lf)
535: {
536: DM_DA *dd = (DM_DA*)da->data;
539: if (lf) *lf = dd->lf;
540: return(0);
541: }
545: /*@C
546: DMDAGetLocalJacobian - Gets from a DM a local jacobian
548: Not Collective
550: Input Parameter:
551: . da - initial distributed array
553: Output Parameter:
554: . lj - the local jacobian
556: Level: intermediate
558: .keywords: distributed array, refine
560: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDAGetLocalFunction(), DMDASetLocalJacobian()
561: @*/
562: PetscErrorCode DMDAGetLocalJacobian(DM da,DMDALocalFunction1 *lj)
563: {
564: DM_DA *dd = (DM_DA*)da->data;
567: if (lj) *lj = dd->lj;
568: return(0);
569: }
573: /*@
574: DMDAComputeFunction - Evaluates a user provided function on each processor that
575: share a DMDA
577: Input Parameters:
578: + da - the DM that defines the grid
579: . vu - input vector
580: . vfu - output vector
581: - w - any user data
583: Notes: Does NOT do ghost updates on vu upon entry
585: This should eventually replace DMDAComputeFunction1
587: Level: advanced
589: .seealso: DMDAComputeJacobian1WithAdic()
591: @*/
592: PetscErrorCode DMDAComputeFunction(DM da,PetscErrorCode (*lf)(void),Vec vu,Vec vfu,void *w)
593: {
595: void *u,*fu;
596: DMDALocalInfo info;
597: PetscErrorCode (*f)(DMDALocalInfo*,void*,void*,void*) = (PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))lf;
598:
600: DMDAGetLocalInfo(da,&info);
601: DMDAVecGetArray(da,vu,&u);
602: DMDAVecGetArray(da,vfu,&fu);
604: (*f)(&info,u,fu,w);
606: DMDAVecRestoreArray(da,vu,&u);
607: DMDAVecRestoreArray(da,vfu,&fu);
608: return(0);
609: }
613: /*@C
614: DMDAComputeFunctionLocal - This is a universal function evaluation routine for
615: a local DM function.
617: Collective on DMDA
619: Input Parameters:
620: + da - the DM context
621: . func - The local function
622: . X - input vector
623: . F - function vector
624: - ctx - A user context
626: Level: intermediate
628: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
629: SNESSetFunction(), SNESSetJacobian()
631: @*/
632: PetscErrorCode DMDAComputeFunctionLocal(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx)
633: {
634: Vec localX;
635: DMDALocalInfo info;
636: void *u;
637: void *fu;
641: DMGetLocalVector(da,&localX);
642: /*
643: Scatter ghost points to local vector, using the 2-step process
644: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
645: */
646: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
647: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
648: DMDAGetLocalInfo(da,&info);
649: DMDAVecGetArray(da,localX,&u);
650: DMDAVecGetArray(da,F,&fu);
651: (*func)(&info,u,fu,ctx);
652: DMDAVecRestoreArray(da,localX,&u);
653: DMDAVecRestoreArray(da,F,&fu);
654: DMRestoreLocalVector(da,&localX);
655: return(0);
656: }
660: /*@C
661: DMDAComputeFunctionLocalGhost - This is a universal function evaluation routine for
662: a local DM function, but the ghost values of the output are communicated and added.
664: Collective on DMDA
666: Input Parameters:
667: + da - the DM context
668: . func - The local function
669: . X - input vector
670: . F - function vector
671: - ctx - A user context
673: Level: intermediate
675: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
676: SNESSetFunction(), SNESSetJacobian()
678: @*/
679: PetscErrorCode DMDAComputeFunctionLocalGhost(DM da, DMDALocalFunction1 func, Vec X, Vec F, void *ctx)
680: {
681: Vec localX, localF;
682: DMDALocalInfo info;
683: void *u;
684: void *fu;
688: DMGetLocalVector(da,&localX);
689: DMGetLocalVector(da,&localF);
690: /*
691: Scatter ghost points to local vector, using the 2-step process
692: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
693: */
694: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
695: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
696: VecSet(F, 0.0);
697: VecSet(localF, 0.0);
698: DMDAGetLocalInfo(da,&info);
699: DMDAVecGetArray(da,localX,&u);
700: DMDAVecGetArray(da,localF,&fu);
701: (*func)(&info,u,fu,ctx);
702: DMLocalToGlobalBegin(da,localF,ADD_VALUES,F);
703: DMLocalToGlobalEnd(da,localF,ADD_VALUES,F);
704: DMDAVecRestoreArray(da,localX,&u);
705: DMDAVecRestoreArray(da,localF,&fu);
706: DMRestoreLocalVector(da,&localX);
707: DMRestoreLocalVector(da,&localF);
708: return(0);
709: }
713: /*@
714: DMDAComputeFunction1 - Evaluates a user provided function on each processor that
715: share a DMDA
717: Input Parameters:
718: + da - the DM that defines the grid
719: . vu - input vector
720: . vfu - output vector
721: - w - any user data
723: Notes: Does NOT do ghost updates on vu upon entry
725: Level: advanced
727: .seealso: DMDAComputeJacobian1WithAdic()
729: @*/
730: PetscErrorCode DMDAComputeFunction1(DM da,Vec vu,Vec vfu,void *w)
731: {
733: void *u,*fu;
734: DMDALocalInfo info;
735: DM_DA *dd = (DM_DA*)da->data;
736:
738: if (!dd->lf) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_NULL,"DMDASetLocalFunction() never called");
739: DMDAGetLocalInfo(da,&info);
740: DMDAVecGetArray(da,vu,&u);
741: DMDAVecGetArray(da,vfu,&fu);
743: CHKMEMQ;
744: (*dd->lf)(&info,u,fu,w);
745: CHKMEMQ;
747: DMDAVecRestoreArray(da,vu,&u);
748: DMDAVecRestoreArray(da,vfu,&fu);
749: return(0);
750: }
754: PetscErrorCode DMDAComputeFunctioniTest1(DM da,void *w)
755: {
756: Vec vu,fu,fui;
758: PetscInt i,n;
759: PetscScalar *ui;
760: PetscRandom rnd;
761: PetscReal norm;
764: DMGetLocalVector(da,&vu);
765: PetscRandomCreate(PETSC_COMM_SELF,&rnd);
766: PetscRandomSetFromOptions(rnd);
767: VecSetRandom(vu,rnd);
768: PetscRandomDestroy(&rnd);
770: DMGetGlobalVector(da,&fu);
771: DMGetGlobalVector(da,&fui);
772:
773: DMDAComputeFunction1(da,vu,fu,w);
775: VecGetArray(fui,&ui);
776: VecGetLocalSize(fui,&n);
777: for (i=0; i<n; i++) {
778: DMDAComputeFunctioni1(da,i,vu,ui+i,w);
779: }
780: VecRestoreArray(fui,&ui);
782: VecAXPY(fui,-1.0,fu);
783: VecNorm(fui,NORM_2,&norm);
784: PetscPrintf(((PetscObject)da)->comm,"Norm of difference in vectors %G\n",norm);
785: VecView(fu,0);
786: VecView(fui,0);
788: DMRestoreLocalVector(da,&vu);
789: DMRestoreGlobalVector(da,&fu);
790: DMRestoreGlobalVector(da,&fui);
791: return(0);
792: }
796: /*@
797: DMDAComputeFunctioni1 - Evaluates a user provided point-wise function
799: Input Parameters:
800: + da - the DM that defines the grid
801: . i - the component of the function we wish to compute (must be local)
802: . vu - input vector
803: . vfu - output value
804: - w - any user data
806: Notes: Does NOT do ghost updates on vu upon entry
808: Level: advanced
810: .seealso: DMDAComputeJacobian1WithAdic()
812: @*/
813: PetscErrorCode DMDAComputeFunctioni1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w)
814: {
816: void *u;
817: DMDALocalInfo info;
818: MatStencil stencil;
819: DM_DA *dd = (DM_DA*)da->data;
820:
823: DMDAGetLocalInfo(da,&info);
824: DMDAVecGetArray(da,vu,&u);
826: /* figure out stencil value from i */
827: stencil.c = i % info.dof;
828: stencil.i = (i % (info.xm*info.dof))/info.dof;
829: stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
830: stencil.k = i/(info.xm*info.ym*info.dof);
832: (*dd->lfi)(&info,&stencil,u,vfu,w);
834: DMDAVecRestoreArray(da,vu,&u);
835: return(0);
836: }
840: /*@
841: DMDAComputeFunctionib1 - Evaluates a user provided point-block function
843: Input Parameters:
844: + da - the DM that defines the grid
845: . i - the component of the function we wish to compute (must be local)
846: . vu - input vector
847: . vfu - output value
848: - w - any user data
850: Notes: Does NOT do ghost updates on vu upon entry
852: Level: advanced
854: .seealso: DMDAComputeJacobian1WithAdic()
856: @*/
857: PetscErrorCode DMDAComputeFunctionib1(DM da,PetscInt i,Vec vu,PetscScalar *vfu,void *w)
858: {
860: void *u;
861: DMDALocalInfo info;
862: MatStencil stencil;
863: DM_DA *dd = (DM_DA*)da->data;
864:
866: DMDAGetLocalInfo(da,&info);
867: DMDAVecGetArray(da,vu,&u);
869: /* figure out stencil value from i */
870: stencil.c = i % info.dof;
871: if (stencil.c) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Point-block functions can only be called for the entire block");
872: stencil.i = (i % (info.xm*info.dof))/info.dof;
873: stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
874: stencil.k = i/(info.xm*info.ym*info.dof);
876: (*dd->lfib)(&info,&stencil,u,vfu,w);
878: DMDAVecRestoreArray(da,vu,&u);
879: return(0);
880: }
882: #if defined(new)
885: /*
886: DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
887: function lives on a DMDA
889: y ~= (F(u + ha) - F(u))/h,
890: where F = nonlinear function, as set by SNESSetFunction()
891: u = current iterate
892: h = difference interval
893: */
894: PetscErrorCode DMDAGetDiagonal_MFFD(DM da,Vec U,Vec a)
895: {
896: PetscScalar h,*aa,*ww,v;
897: PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
899: PetscInt gI,nI;
900: MatStencil stencil;
901: DMDALocalInfo info;
902:
904: (*ctx->func)(0,U,a,ctx->funcctx);
905: (*ctx->funcisetbase)(U,ctx->funcctx);
907: VecGetArray(U,&ww);
908: VecGetArray(a,&aa);
909:
910: nI = 0;
911: h = ww[gI];
912: if (h == 0.0) h = 1.0;
913: #if !defined(PETSC_USE_COMPLEX)
914: if (h < umin && h >= 0.0) h = umin;
915: else if (h < 0.0 && h > -umin) h = -umin;
916: #else
917: if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin;
918: else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
919: #endif
920: h *= epsilon;
921:
922: ww[gI] += h;
923: (*ctx->funci)(i,w,&v,ctx->funcctx);
924: aa[nI] = (v - aa[nI])/h;
925: ww[gI] -= h;
926: nI++;
927: }
928: VecRestoreArray(U,&ww);
929: VecRestoreArray(a,&aa);
930: return(0);
931: }
932: #endif
934: #if defined(PETSC_HAVE_ADIC)
935: EXTERN_C_BEGIN
936: #include <adic/ad_utils.h>
937: EXTERN_C_END
941: /*@C
942: DMDAComputeJacobian1WithAdic - Evaluates a adiC provided Jacobian function on each processor that
943: share a DMDA
945: Input Parameters:
946: + da - the DM that defines the grid
947: . vu - input vector (ghosted)
948: . J - output matrix
949: - w - any user data
951: Level: advanced
953: Notes: Does NOT do ghost updates on vu upon entry
955: .seealso: DMDAComputeFunction1()
957: @*/
958: PetscErrorCode DMDAComputeJacobian1WithAdic(DM da,Vec vu,Mat J,void *w)
959: {
961: PetscInt gtdof,tdof;
962: PetscScalar *ustart;
963: DMDALocalInfo info;
964: void *ad_u,*ad_f,*ad_ustart,*ad_fstart;
965: ISColoring iscoloring;
968: DMDAGetLocalInfo(da,&info);
970: PetscADResetIndep();
972: /* get space for derivative objects. */
973: DMDAGetAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,>dof);
974: DMDAGetAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
975: VecGetArray(vu,&ustart);
976: DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
978: PetscADSetValueAndColor(ad_ustart,gtdof,iscoloring->colors,ustart);
980: VecRestoreArray(vu,&ustart);
981: ISColoringDestroy(&iscoloring);
982: PetscADIncrementTotalGradSize(iscoloring->n);
983: PetscADSetIndepDone();
985: PetscLogEventBegin(DMDA_LocalADFunction,0,0,0,0);
986: (*dd->adic_lf)(&info,ad_u,ad_f,w);
987: PetscLogEventEnd(DMDA_LocalADFunction,0,0,0,0);
989: /* stick the values into the matrix */
990: MatSetValuesAdic(J,(PetscScalar**)ad_fstart);
992: /* return space for derivative objects. */
993: DMDARestoreAdicArray(da,PETSC_TRUE,&ad_u,&ad_ustart,>dof);
994: DMDARestoreAdicArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
995: return(0);
996: }
1000: /*@C
1001: DMDAMultiplyByJacobian1WithAdic - Applies an ADIC-provided Jacobian function to a vector on
1002: each processor that shares a DMDA.
1004: Input Parameters:
1005: + da - the DM that defines the grid
1006: . vu - Jacobian is computed at this point (ghosted)
1007: . v - product is done on this vector (ghosted)
1008: . fu - output vector = J(vu)*v (not ghosted)
1009: - w - any user data
1011: Notes:
1012: This routine does NOT do ghost updates on vu upon entry.
1014: Level: advanced
1016: .seealso: DMDAComputeFunction1()
1018: @*/
1019: PetscErrorCode DMDAMultiplyByJacobian1WithAdic(DM da,Vec vu,Vec v,Vec f,void *w)
1020: {
1022: PetscInt i,gtdof,tdof;
1023: PetscScalar *avu,*av,*af,*ad_vustart,*ad_fstart;
1024: DMDALocalInfo info;
1025: void *ad_vu,*ad_f;
1028: DMDAGetLocalInfo(da,&info);
1030: /* get space for derivative objects. */
1031: DMDAGetAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,>dof);
1032: DMDAGetAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
1034: /* copy input vector into derivative object */
1035: VecGetArray(vu,&avu);
1036: VecGetArray(v,&av);
1037: for (i=0; i<gtdof; i++) {
1038: ad_vustart[2*i] = avu[i];
1039: ad_vustart[2*i+1] = av[i];
1040: }
1041: VecRestoreArray(vu,&avu);
1042: VecRestoreArray(v,&av);
1044: PetscADResetIndep();
1045: PetscADIncrementTotalGradSize(1);
1046: PetscADSetIndepDone();
1048: (*dd->adicmf_lf)(&info,ad_vu,ad_f,w);
1050: /* stick the values into the vector */
1051: VecGetArray(f,&af);
1052: for (i=0; i<tdof; i++) {
1053: af[i] = ad_fstart[2*i+1];
1054: }
1055: VecRestoreArray(f,&af);
1057: /* return space for derivative objects. */
1058: DMDARestoreAdicMFArray(da,PETSC_TRUE,&ad_vu,&ad_vustart,>dof);
1059: DMDARestoreAdicMFArray(da,PETSC_FALSE,&ad_f,&ad_fstart,&tdof);
1060: return(0);
1061: }
1062: #endif
1066: /*@
1067: DMDAComputeJacobian1 - Evaluates a local Jacobian function on each processor that
1068: share a DMDA
1070: Input Parameters:
1071: + da - the DM that defines the grid
1072: . vu - input vector (ghosted)
1073: . J - output matrix
1074: - w - any user data
1076: Notes: Does NOT do ghost updates on vu upon entry
1078: Level: advanced
1080: .seealso: DMDAComputeFunction1()
1082: @*/
1083: PetscErrorCode DMDAComputeJacobian1(DM da,Vec vu,Mat J,void *w)
1084: {
1086: void *u;
1087: DMDALocalInfo info;
1088: DM_DA *dd = (DM_DA*)da->data;
1091: DMDAGetLocalInfo(da,&info);
1092: DMDAVecGetArray(da,vu,&u);
1093: (*dd->lj)(&info,u,J,w);
1094: DMDAVecRestoreArray(da,vu,&u);
1095: return(0);
1096: }
1101: /*
1102: DMDAComputeJacobian1WithAdifor - Evaluates a ADIFOR provided Jacobian local function on each processor that
1103: share a DMDA
1105: Input Parameters:
1106: + da - the DM that defines the grid
1107: . vu - input vector (ghosted)
1108: . J - output matrix
1109: - w - any user data
1111: Notes: Does NOT do ghost updates on vu upon entry
1113: .seealso: DMDAComputeFunction1()
1115: */
1116: PetscErrorCode DMDAComputeJacobian1WithAdifor(DM da,Vec vu,Mat J,void *w)
1117: {
1118: PetscErrorCode ierr;
1119: PetscInt i,Nc,N;
1120: ISColoringValue *color;
1121: DMDALocalInfo info;
1122: PetscScalar *u,*g_u,*g_f,*f = 0,*p_u;
1123: ISColoring iscoloring;
1124: DM_DA *dd = (DM_DA*)da->data;
1125: void (*lf)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*) =
1126: (void (*)(PetscInt*,DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*))*dd->adifor_lf;
1129: DMCreateColoring(da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
1130: Nc = iscoloring->n;
1131: DMDAGetLocalInfo(da,&info);
1132: N = info.gxm*info.gym*info.gzm*info.dof;
1134: /* get space for derivative objects. */
1135: PetscMalloc(Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar),&g_u);
1136: PetscMemzero(g_u,Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar));
1137: p_u = g_u;
1138: color = iscoloring->colors;
1139: for (i=0; i<N; i++) {
1140: p_u[*color++] = 1.0;
1141: p_u += Nc;
1142: }
1143: ISColoringDestroy(&iscoloring);
1144: PetscMalloc2(Nc*info.xm*info.ym*info.zm*info.dof,PetscScalar,&g_f,info.xm*info.ym*info.zm*info.dof,PetscScalar,&f);
1146: /* Seed the input array g_u with coloring information */
1147:
1148: VecGetArray(vu,&u);
1149: (lf)(&Nc,&info,u,g_u,&Nc,f,g_f,&Nc,w,&ierr);
1150: VecRestoreArray(vu,&u);
1152: /* stick the values into the matrix */
1153: /* PetscScalarView(Nc*info.xm*info.ym,g_f,0); */
1154: MatSetValuesAdifor(J,Nc,g_f);
1156: /* return space for derivative objects. */
1157: PetscFree(g_u);
1158: PetscFree2(g_f,f);
1159: return(0);
1160: }
1164: /*@C
1165: DMDAFormjacobianLocal - This is a universal Jacobian evaluation routine for
1166: a local DM function.
1168: Collective on DMDA
1170: Input Parameters:
1171: + da - the DM context
1172: . func - The local function
1173: . X - input vector
1174: . J - Jacobian matrix
1175: - ctx - A user context
1177: Level: intermediate
1179: .seealso: DMDASetLocalFunction(), DMDASetLocalJacobian(), DMDASetLocalAdicFunction(), DMDASetLocalAdicMFFunction(),
1180: SNESSetFunction(), SNESSetJacobian()
1182: @*/
1183: PetscErrorCode DMDAFormJacobianLocal(DM da, DMDALocalFunction1 func, Vec X, Mat J, void *ctx)
1184: {
1185: Vec localX;
1186: DMDALocalInfo info;
1187: void *u;
1191: DMGetLocalVector(da,&localX);
1192: /*
1193: Scatter ghost points to local vector, using the 2-step process
1194: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
1195: */
1196: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
1197: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
1198: DMDAGetLocalInfo(da,&info);
1199: DMDAVecGetArray(da,localX,&u);
1200: (*func)(&info,u,J,ctx);
1201: DMDAVecRestoreArray(da,localX,&u);
1202: DMRestoreLocalVector(da,&localX);
1203: return(0);
1204: }
1208: /*@C
1209: DMDAMultiplyByJacobian1WithAD - Applies a Jacobian function supplied by ADIFOR or ADIC
1210: to a vector on each processor that shares a DMDA.
1212: Input Parameters:
1213: + da - the DM that defines the grid
1214: . vu - Jacobian is computed at this point (ghosted)
1215: . v - product is done on this vector (ghosted)
1216: . fu - output vector = J(vu)*v (not ghosted)
1217: - w - any user data
1219: Notes:
1220: This routine does NOT do ghost updates on vu and v upon entry.
1221:
1222: Automatically calls DMDAMultiplyByJacobian1WithAdifor() or DMDAMultiplyByJacobian1WithAdic()
1223: depending on whether DMDASetLocalAdicMFFunction() or DMDASetLocalAdiforMFFunction() was called.
1225: Level: advanced
1227: .seealso: DMDAComputeFunction1(), DMDAMultiplyByJacobian1WithAdifor(), DMDAMultiplyByJacobian1WithAdic()
1229: @*/
1230: PetscErrorCode DMDAMultiplyByJacobian1WithAD(DM da,Vec u,Vec v,Vec f,void *w)
1231: {
1233: DM_DA *dd = (DM_DA*)da->data;
1236: if (dd->adicmf_lf) {
1237: #if defined(PETSC_HAVE_ADIC)
1238: DMDAMultiplyByJacobian1WithAdic(da,u,v,f,w);
1239: #else
1240: SETERRQ(((PetscObject)da)->comm,PETSC_ERR_SUP_SYS,"Requires ADIC to be installed and cannot use complex numbers");
1241: #endif
1242: } else if (dd->adiformf_lf) {
1243: DMDAMultiplyByJacobian1WithAdifor(da,u,v,f,w);
1244: } else SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ORDER,"Must call DMDASetLocalAdiforMFFunction() or DMDASetLocalAdicMFFunction() before using");
1245: return(0);
1246: }
1251: /*@C
1252: DMDAMultiplyByJacobian1WithAdifor - Applies a ADIFOR provided Jacobian function on each processor that
1253: share a DM to a vector
1255: Input Parameters:
1256: + da - the DM that defines the grid
1257: . vu - Jacobian is computed at this point (ghosted)
1258: . v - product is done on this vector (ghosted)
1259: . fu - output vector = J(vu)*v (not ghosted)
1260: - w - any user data
1262: Notes: Does NOT do ghost updates on vu and v upon entry
1264: Level: advanced
1266: .seealso: DMDAComputeFunction1()
1268: @*/
1269: PetscErrorCode DMDAMultiplyByJacobian1WithAdifor(DM da,Vec u,Vec v,Vec f,void *w)
1270: {
1272: PetscScalar *au,*av,*af,*awork;
1273: Vec work;
1274: DMDALocalInfo info;
1275: DM_DA *dd = (DM_DA*)da->data;
1276: void (*lf)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*) =
1277: (void (*)(DMDALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*))*dd->adiformf_lf;
1280: DMDAGetLocalInfo(da,&info);
1282: DMGetGlobalVector(da,&work);
1283: VecGetArray(u,&au);
1284: VecGetArray(v,&av);
1285: VecGetArray(f,&af);
1286: VecGetArray(work,&awork);
1287: (lf)(&info,au,av,awork,af,w,&ierr);
1288: VecRestoreArray(u,&au);
1289: VecRestoreArray(v,&av);
1290: VecRestoreArray(f,&af);
1291: VecRestoreArray(work,&awork);
1292: DMRestoreGlobalVector(da,&work);
1294: return(0);
1295: }
1299: PetscErrorCode DMSetUp_DA_2D(DM da)
1300: {
1301: DM_DA *dd = (DM_DA*)da->data;
1302: const PetscInt M = dd->M;
1303: const PetscInt N = dd->N;
1304: PetscInt m = dd->m;
1305: PetscInt n = dd->n;
1306: const PetscInt dof = dd->w;
1307: const PetscInt s = dd->s;
1308: const DMDABoundaryType bx = dd->bx;
1309: const DMDABoundaryType by = dd->by;
1310: const DMDAStencilType stencil_type = dd->stencil_type;
1311: PetscInt *lx = dd->lx;
1312: PetscInt *ly = dd->ly;
1313: MPI_Comm comm;
1314: PetscMPIInt rank,size;
1315: PetscInt xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end,IXs,IXe,IYs,IYe;
1316: PetscInt up,down,left,right,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn,*idx_cpy;
1317: const PetscInt *idx_full;
1318: PetscInt xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
1319: PetscInt s_x,s_y; /* s proportionalized to w */
1320: PetscInt sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0;
1321: Vec local,global;
1322: VecScatter ltog,gtol;
1323: IS to,from,ltogis;
1324: PetscErrorCode ierr;
1327: if (dof < 1) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
1328: if (s < 0) SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
1329: PetscObjectGetComm((PetscObject)da,&comm);
1330: #if !defined(PETSC_USE_64BIT_INDICES)
1331: if (((Petsc64bitInt) M)*((Petsc64bitInt) N)*((Petsc64bitInt) dof) > (Petsc64bitInt) PETSC_MPI_INT_MAX) SETERRQ3(comm,PETSC_ERR_INT_OVERFLOW,"Mesh of %D by %D by %D (dof) is too large for 32 bit indices",M,N,dof);
1332: #endif
1334: if (dof < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
1335: if (s < 0) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);
1337: MPI_Comm_size(comm,&size);
1338: MPI_Comm_rank(comm,&rank);
1340: dd->dim = 2;
1341: PetscMalloc(dof*sizeof(char*),&dd->fieldname);
1342: PetscMemzero(dd->fieldname,dof*sizeof(char*));
1344: if (m != PETSC_DECIDE) {
1345: if (m < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);
1346: else if (m > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);
1347: }
1348: if (n != PETSC_DECIDE) {
1349: if (n < 1) SETERRQ1(comm,PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);
1350: else if (n > size) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);
1351: }
1353: if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
1354: if (n != PETSC_DECIDE) {
1355: m = size/n;
1356: } else if (m != PETSC_DECIDE) {
1357: n = size/m;
1358: } else {
1359: /* try for squarish distribution */
1360: m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N)));
1361: if (!m) m = 1;
1362: while (m > 0) {
1363: n = size/m;
1364: if (m*n == size) break;
1365: m--;
1366: }
1367: if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
1368: }
1369: if (m*n != size) SETERRQ(comm,PETSC_ERR_PLIB,"Unable to create partition, check the size of the communicator and input m and n ");
1370: } else if (m*n != size) SETERRQ(comm,PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");
1372: if (M < m) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
1373: if (N < n) SETERRQ2(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);
1375: /*
1376: Determine locally owned region
1377: xs is the first local node number, x is the number of local nodes
1378: */
1379: if (!lx) {
1380: PetscMalloc(m*sizeof(PetscInt), &dd->lx);
1381: lx = dd->lx;
1382: for (i=0; i<m; i++) {
1383: lx[i] = M/m + ((M % m) > i);
1384: }
1385: }
1386: x = lx[rank % m];
1387: xs = 0;
1388: for (i=0; i<(rank % m); i++) {
1389: xs += lx[i];
1390: }
1391: #if defined(PETSC_USE_DEBUG)
1392: left = xs;
1393: for (i=(rank % m); i<m; i++) {
1394: left += lx[i];
1395: }
1396: if (left != M) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
1397: #endif
1399: /*
1400: Determine locally owned region
1401: ys is the first local node number, y is the number of local nodes
1402: */
1403: if (!ly) {
1404: PetscMalloc(n*sizeof(PetscInt), &dd->ly);
1405: ly = dd->ly;
1406: for (i=0; i<n; i++) {
1407: ly[i] = N/n + ((N % n) > i);
1408: }
1409: }
1410: y = ly[rank/m];
1411: ys = 0;
1412: for (i=0; i<(rank/m); i++) {
1413: ys += ly[i];
1414: }
1415: #if defined(PETSC_USE_DEBUG)
1416: left = ys;
1417: for (i=(rank/m); i<n; i++) {
1418: left += ly[i];
1419: }
1420: if (left != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
1421: #endif
1423: /*
1424: check if the scatter requires more than one process neighbor or wraps around
1425: the domain more than once
1426: */
1427: if ((x < s) && ((m > 1) || (bx == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local x-width of domain x %D is smaller than stencil width s %D",x,s);
1428: if ((y < s) && ((n > 1) || (by == DMDA_BOUNDARY_PERIODIC))) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local y-width of domain y %D is smaller than stencil width s %D",y,s);
1429: xe = xs + x;
1430: ye = ys + y;
1432: /* determine ghost region (Xs) and region scattered into (IXs) */
1433: /* Assume No Periodicity */
1434: if (xs-s > 0) { Xs = xs - s; IXs = xs - s; } else { Xs = 0; IXs = 0; }
1435: if (xe+s <= M) { Xe = xe + s; IXe = xe + s; } else { Xe = M; IXe = M; }
1436: if (ys-s > 0) { Ys = ys - s; IYs = ys - s; } else { Ys = 0; IYs = 0; }
1437: if (ye+s <= N) { Ye = ye + s; IYe = ye + s; } else { Ye = N; IYe = N; }
1439: /* fix for periodicity/ghosted */
1440: if (bx) { Xs = xs - s; Xe = xe + s; }
1441: if (bx == DMDA_BOUNDARY_PERIODIC) { IXs = xs - s; IXe = xe + s; }
1442: if (by) { Ys = ys - s; Ye = ye + s; }
1443: if (by == DMDA_BOUNDARY_PERIODIC) { IYs = ys - s; IYe = ye + s; }
1445: /* Resize all X parameters to reflect w */
1446: s_x = s;
1447: s_y = s;
1449: /* determine starting point of each processor */
1450: nn = x*y;
1451: PetscMalloc2(size+1,PetscInt,&bases,size,PetscInt,&ldims);
1452: MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
1453: bases[0] = 0;
1454: for (i=1; i<=size; i++) {
1455: bases[i] = ldims[i-1];
1456: }
1457: for (i=1; i<=size; i++) {
1458: bases[i] += bases[i-1];
1459: }
1460: base = bases[rank]*dof;
1462: /* allocate the base parallel and sequential vectors */
1463: dd->Nlocal = x*y*dof;
1464: VecCreateMPIWithArray(comm,dof,dd->Nlocal,PETSC_DECIDE,0,&global);
1465: dd->nlocal = (Xe-Xs)*(Ye-Ys)*dof;
1466: VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->nlocal,0,&local);
1468: /* generate appropriate vector scatters */
1469: /* local to global inserts non-ghost point region into global */
1470: VecGetOwnershipRange(global,&start,&end);
1471: ISCreateStride(comm,x*y*dof,start,1,&to);
1473: count = x*y;
1474: PetscMalloc(x*y*sizeof(PetscInt),&idx);
1475: left = xs - Xs; right = left + x;
1476: down = ys - Ys; up = down + y;
1477: count = 0;
1478: for (i=down; i<up; i++) {
1479: for (j=left; j<right; j++) {
1480: idx[count++] = i*(Xe-Xs) + j;
1481: }
1482: }
1484: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&from);
1485: VecScatterCreate(local,from,global,to,<og);
1486: PetscLogObjectParent(dd,ltog);
1487: ISDestroy(&from);
1488: ISDestroy(&to);
1490: /* global to local must include ghost points within the domain,
1491: but not ghost points outside the domain that aren't periodic */
1492: if (stencil_type == DMDA_STENCIL_BOX) {
1493: count = (IXe-IXs)*(IYe-IYs);
1494: PetscMalloc(count*sizeof(PetscInt),&idx);
1496: left = IXs - Xs; right = left + (IXe-IXs);
1497: down = IYs - Ys; up = down + (IYe-IYs);
1498: count = 0;
1499: for (i=down; i<up; i++) {
1500: for (j=left; j<right; j++) {
1501: idx[count++] = j + i*(Xe-Xs);
1502: }
1503: }
1504: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
1506: } else {
1507: /* must drop into cross shape region */
1508: /* ---------|
1509: | top |
1510: |--- ---| up
1511: | middle |
1512: | |
1513: ---- ---- down
1514: | bottom |
1515: -----------
1516: Xs xs xe Xe */
1517: count = (ys-IYs)*x + y*(IXe-IXs) + (IYe-ye)*x;
1518: PetscMalloc(count*sizeof(PetscInt),&idx);
1520: left = xs - Xs; right = left + x;
1521: down = ys - Ys; up = down + y;
1522: count = 0;
1523: /* bottom */
1524: for (i=(IYs-Ys); i<down; i++) {
1525: for (j=left; j<right; j++) {
1526: idx[count++] = j + i*(Xe-Xs);
1527: }
1528: }
1529: /* middle */
1530: for (i=down; i<up; i++) {
1531: for (j=(IXs-Xs); j<(IXe-Xs); j++) {
1532: idx[count++] = j + i*(Xe-Xs);
1533: }
1534: }
1535: /* top */
1536: for (i=up; i<up+IYe-ye; i++) {
1537: for (j=left; j<right; j++) {
1538: idx[count++] = j + i*(Xe-Xs);
1539: }
1540: }
1541: ISCreateBlock(comm,dof,count,idx,PETSC_OWN_POINTER,&to);
1542: }
1545: /* determine who lies on each side of us stored in n6 n7 n8
1546: n3 n5
1547: n0 n1 n2
1548: */
1550: /* Assume the Non-Periodic Case */
1551: n1 = rank - m;
1552: if (rank % m) {
1553: n0 = n1 - 1;
1554: } else {
1555: n0 = -1;
1556: }
1557: if ((rank+1) % m) {
1558: n2 = n1 + 1;
1559: n5 = rank + 1;
1560: n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
1561: } else {
1562: n2 = -1; n5 = -1; n8 = -1;
1563: }
1564: if (rank % m) {
1565: n3 = rank - 1;
1566: n6 = n3 + m; if (n6 >= m*n) n6 = -1;
1567: } else {
1568: n3 = -1; n6 = -1;
1569: }
1570: n7 = rank + m; if (n7 >= m*n) n7 = -1;
1572: if (bx == DMDA_BOUNDARY_PERIODIC && by == DMDA_BOUNDARY_PERIODIC) {
1573: /* Modify for Periodic Cases */
1574: /* Handle all four corners */
1575: if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
1576: if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
1577: if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
1578: if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;
1580: /* Handle Top and Bottom Sides */
1581: if (n1 < 0) n1 = rank + m * (n-1);
1582: if (n7 < 0) n7 = rank - m * (n-1);
1583: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
1584: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
1585: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
1586: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
1588: /* Handle Left and Right Sides */
1589: if (n3 < 0) n3 = rank + (m-1);
1590: if (n5 < 0) n5 = rank - (m-1);
1591: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
1592: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
1593: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
1594: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
1595: } else if (by == DMDA_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */
1596: if (n1 < 0) n1 = rank + m * (n-1);
1597: if (n7 < 0) n7 = rank - m * (n-1);
1598: if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
1599: if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
1600: if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
1601: if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
1602: } else if (bx == DMDA_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */
1603: if (n3 < 0) n3 = rank + (m-1);
1604: if (n5 < 0) n5 = rank - (m-1);
1605: if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
1606: if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
1607: if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
1608: if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
1609: }
1611: PetscMalloc(9*sizeof(PetscInt),&dd->neighbors);
1612: dd->neighbors[0] = n0;
1613: dd->neighbors[1] = n1;
1614: dd->neighbors[2] = n2;
1615: dd->neighbors[3] = n3;
1616: dd->neighbors[4] = rank;
1617: dd->neighbors[5] = n5;
1618: dd->neighbors[6] = n6;
1619: dd->neighbors[7] = n7;
1620: dd->neighbors[8] = n8;
1622: if (stencil_type == DMDA_STENCIL_STAR) {
1623: /* save corner processor numbers */
1624: sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
1625: n0 = n2 = n6 = n8 = -1;
1626: }
1628: PetscMalloc((Xe-Xs)*(Ye-Ys)*sizeof(PetscInt),&idx);
1629: PetscLogObjectMemory(da,(Xe-Xs)*(Ye-Ys)*sizeof(PetscInt));
1631: nn = 0;
1632: xbase = bases[rank];
1633: for (i=1; i<=s_y; i++) {
1634: if (n0 >= 0) { /* left below */
1635: x_t = lx[n0 % m];
1636: y_t = ly[(n0/m)];
1637: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
1638: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1639: }
1640: if (n1 >= 0) { /* directly below */
1641: x_t = x;
1642: y_t = ly[(n1/m)];
1643: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
1644: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1645: }
1646: if (n2 >= 0) { /* right below */
1647: x_t = lx[n2 % m];
1648: y_t = ly[(n2/m)];
1649: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
1650: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1651: }
1652: }
1654: for (i=0; i<y; i++) {
1655: if (n3 >= 0) { /* directly left */
1656: x_t = lx[n3 % m];
1657: /* y_t = y; */
1658: s_t = bases[n3] + (i+1)*x_t - s_x;
1659: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1660: }
1662: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
1664: if (n5 >= 0) { /* directly right */
1665: x_t = lx[n5 % m];
1666: /* y_t = y; */
1667: s_t = bases[n5] + (i)*x_t;
1668: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1669: }
1670: }
1672: for (i=1; i<=s_y; i++) {
1673: if (n6 >= 0) { /* left above */
1674: x_t = lx[n6 % m];
1675: /* y_t = ly[(n6/m)]; */
1676: s_t = bases[n6] + (i)*x_t - s_x;
1677: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1678: }
1679: if (n7 >= 0) { /* directly above */
1680: x_t = x;
1681: /* y_t = ly[(n7/m)]; */
1682: s_t = bases[n7] + (i-1)*x_t;
1683: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1684: }
1685: if (n8 >= 0) { /* right above */
1686: x_t = lx[n8 % m];
1687: /* y_t = ly[(n8/m)]; */
1688: s_t = bases[n8] + (i-1)*x_t;
1689: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1690: }
1691: }
1693: ISCreateBlock(comm,dof,nn,idx,PETSC_COPY_VALUES,&from);
1694: VecScatterCreate(global,from,local,to,>ol);
1695: PetscLogObjectParent(da,gtol);
1696: ISDestroy(&to);
1697: ISDestroy(&from);
1699: if (stencil_type == DMDA_STENCIL_STAR) {
1700: n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
1701: }
1703: if ((stencil_type == DMDA_STENCIL_STAR) ||
1704: (bx && bx != DMDA_BOUNDARY_PERIODIC) ||
1705: (by && by != DMDA_BOUNDARY_PERIODIC)) {
1706: /*
1707: Recompute the local to global mappings, this time keeping the
1708: information about the cross corner processor numbers and any ghosted
1709: but not periodic indices.
1710: */
1711: nn = 0;
1712: xbase = bases[rank];
1713: for (i=1; i<=s_y; i++) {
1714: if (n0 >= 0) { /* left below */
1715: x_t = lx[n0 % m];
1716: y_t = ly[(n0/m)];
1717: s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
1718: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1719: } else if (xs-Xs > 0 && ys-Ys > 0) {
1720: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1721: }
1722: if (n1 >= 0) { /* directly below */
1723: x_t = x;
1724: y_t = ly[(n1/m)];
1725: s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
1726: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1727: } else if (ys-Ys > 0) {
1728: for (j=0; j<x; j++) { idx[nn++] = -1;}
1729: }
1730: if (n2 >= 0) { /* right below */
1731: x_t = lx[n2 % m];
1732: y_t = ly[(n2/m)];
1733: s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
1734: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1735: } else if (Xe-xe> 0 && ys-Ys > 0) {
1736: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1737: }
1738: }
1740: for (i=0; i<y; i++) {
1741: if (n3 >= 0) { /* directly left */
1742: x_t = lx[n3 % m];
1743: /* y_t = y; */
1744: s_t = bases[n3] + (i+1)*x_t - s_x;
1745: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1746: } else if (xs-Xs > 0) {
1747: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1748: }
1750: for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */
1752: if (n5 >= 0) { /* directly right */
1753: x_t = lx[n5 % m];
1754: /* y_t = y; */
1755: s_t = bases[n5] + (i)*x_t;
1756: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1757: } else if (Xe-xe > 0) {
1758: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1759: }
1760: }
1762: for (i=1; i<=s_y; i++) {
1763: if (n6 >= 0) { /* left above */
1764: x_t = lx[n6 % m];
1765: /* y_t = ly[(n6/m)]; */
1766: s_t = bases[n6] + (i)*x_t - s_x;
1767: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1768: } else if (xs-Xs > 0 && Ye-ye > 0) {
1769: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1770: }
1771: if (n7 >= 0) { /* directly above */
1772: x_t = x;
1773: /* y_t = ly[(n7/m)]; */
1774: s_t = bases[n7] + (i-1)*x_t;
1775: for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
1776: } else if (Ye-ye > 0) {
1777: for (j=0; j<x; j++) { idx[nn++] = -1;}
1778: }
1779: if (n8 >= 0) { /* right above */
1780: x_t = lx[n8 % m];
1781: /* y_t = ly[(n8/m)]; */
1782: s_t = bases[n8] + (i-1)*x_t;
1783: for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
1784: } else if (Xe-xe > 0 && Ye-ye > 0) {
1785: for (j=0; j<s_x; j++) { idx[nn++] = -1;}
1786: }
1787: }
1788: }
1789: /*
1790: Set the local to global ordering in the global vector, this allows use
1791: of VecSetValuesLocal().
1792: */
1793: ISCreateBlock(comm,dof,nn,idx,PETSC_OWN_POINTER,<ogis);
1794: PetscMalloc(nn*dof*sizeof(PetscInt),&idx_cpy);
1795: PetscLogObjectMemory(da,nn*dof*sizeof(PetscInt));
1796: ISGetIndices(ltogis, &idx_full);
1797: PetscMemcpy(idx_cpy,idx_full,nn*dof*sizeof(PetscInt));
1798: ISRestoreIndices(ltogis, &idx_full);
1799: ISLocalToGlobalMappingCreateIS(ltogis,&da->ltogmap);
1800: PetscLogObjectParent(da,da->ltogmap);
1801: ISDestroy(<ogis);
1802: ISLocalToGlobalMappingBlock(da->ltogmap,dd->w,&da->ltogmapb);
1803: PetscLogObjectParent(da,da->ltogmap);
1805: PetscFree2(bases,ldims);
1806: dd->m = m; dd->n = n;
1807: /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */
1808: dd->xs = xs*dof; dd->xe = xe*dof; dd->ys = ys; dd->ye = ye; dd->zs = 0; dd->ze = 1;
1809: dd->Xs = Xs*dof; dd->Xe = Xe*dof; dd->Ys = Ys; dd->Ye = Ye; dd->Zs = 0; dd->Ze = 1;
1811: VecDestroy(&local);
1812: VecDestroy(&global);
1814: dd->gtol = gtol;
1815: dd->ltog = ltog;
1816: dd->idx = idx_cpy;
1817: dd->Nl = nn*dof;
1818: dd->base = base;
1819: da->ops->view = DMView_DA_2d;
1820: dd->ltol = PETSC_NULL;
1821: dd->ao = PETSC_NULL;
1823: return(0);
1824: }
1828: /*@C
1829: DMDACreate2d - Creates an object that will manage the communication of two-dimensional
1830: regular array data that is distributed across some processors.
1832: Collective on MPI_Comm
1834: Input Parameters:
1835: + comm - MPI communicator
1836: . bx,by - type of ghost nodes the array have.
1837: Use one of DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_GHOSTED, DMDA_BOUNDARY_PERIODIC.
1838: . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
1839: . M,N - global dimension in each direction of the array (use -M and or -N to indicate that it may be set to a different value
1840: from the command line with -da_grid_x <M> -da_grid_y <N>)
1841: . m,n - corresponding number of processors in each dimension
1842: (or PETSC_DECIDE to have calculated)
1843: . dof - number of degrees of freedom per node
1844: . s - stencil width
1845: - lx, ly - arrays containing the number of nodes in each cell along
1846: the x and y coordinates, or PETSC_NULL. If non-null, these
1847: must be of length as m and n, and the corresponding
1848: m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
1849: must be M, and the sum of the ly[] entries must be N.
1851: Output Parameter:
1852: . da - the resulting distributed array object
1854: Options Database Key:
1855: + -da_view - Calls DMView() at the conclusion of DMDACreate2d()
1856: . -da_grid_x <nx> - number of grid points in x direction, if M < 0
1857: . -da_grid_y <ny> - number of grid points in y direction, if N < 0
1858: . -da_processors_x <nx> - number of processors in x direction
1859: . -da_processors_y <ny> - number of processors in y direction
1860: . -da_refine_x <rx> - refinement ratio in x direction
1861: . -da_refine_y <ry> - refinement ratio in y direction
1862: - -da_refine <n> - refine the DMDA n times before creating, if M or N < 0
1865: Level: beginner
1867: Notes:
1868: The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the
1869: standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes
1870: the standard 9-pt stencil.
1872: The array data itself is NOT stored in the DMDA, it is stored in Vec objects;
1873: The appropriate vector objects can be obtained with calls to DMCreateGlobalVector()
1874: and DMCreateLocalVector() and calls to VecDuplicate() if more are needed.
1876: .keywords: distributed array, create, two-dimensional
1878: .seealso: DMDestroy(), DMView(), DMDACreate1d(), DMDACreate3d(), DMGlobalToLocalBegin(), DMDAGetRefinementFactor(),
1879: DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMDALocalToLocalBegin(), DMDALocalToLocalEnd(), DMDASetRefinementFactor(),
1880: DMDAGetInfo(), DMCreateGlobalVector(), DMCreateLocalVector(), DMDACreateNaturalVector(), DMLoad(), DMDAGetOwnershipRanges()
1882: @*/
1884: PetscErrorCode DMDACreate2d(MPI_Comm comm,DMDABoundaryType bx,DMDABoundaryType by,DMDAStencilType stencil_type,
1885: PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,const PetscInt lx[],const PetscInt ly[],DM *da)
1886: {
1890: DMDACreate(comm, da);
1891: DMDASetDim(*da, 2);
1892: DMDASetSizes(*da, M, N, 1);
1893: DMDASetNumProcs(*da, m, n, PETSC_DECIDE);
1894: DMDASetBoundaryType(*da, bx, by, DMDA_BOUNDARY_NONE);
1895: DMDASetDof(*da, dof);
1896: DMDASetStencilType(*da, stencil_type);
1897: DMDASetStencilWidth(*da, s);
1898: DMDASetOwnershipRanges(*da, lx, ly, PETSC_NULL);
1899: /* This violates the behavior for other classes, but right now users expect negative dimensions to be handled this way */
1900: DMSetFromOptions(*da);
1901: DMSetUp(*da);
1902: DMView_DA_Private(*da);
1903: return(0);
1904: }