Actual source code: wb.c
2: #include <petscdmda.h>
3: #include <petsc/private/pcmgimpl.h>
4: #include <petscctable.h>
6: typedef struct {
7: PCExoticType type;
8: Mat P; /* the constructed interpolation matrix */
9: PetscBool directSolve; /* use direct LU factorization to construct interpolation */
10: KSP ksp;
11: } PC_Exotic;
13: const char *const PCExoticTypes[] = {"face","wirebasket","PCExoticType","PC_Exotic",NULL};
15: /*
16: DMDAGetWireBasketInterpolation - Gets the interpolation for a wirebasket based coarse space
18: */
19: PetscErrorCode DMDAGetWireBasketInterpolation(PC pc,DM da,PC_Exotic *exotic,Mat Aglobal,MatReuse reuse,Mat *P)
20: {
21: PetscInt dim,i,j,k,m,n,p,dof,Nint,Nface,Nwire,Nsurf,*Iint,*Isurf,cint = 0,csurf = 0,istart,jstart,kstart,*II,N,c = 0;
22: PetscInt mwidth,nwidth,pwidth,cnt,mp,np,pp,Ntotal,gl[26],*globals,Ng,*IIint,*IIsurf,Nt;
23: Mat Xint, Xsurf,Xint_tmp;
24: IS isint,issurf,is,row,col;
25: ISLocalToGlobalMapping ltg;
26: MPI_Comm comm;
27: Mat A,Aii,Ais,Asi,*Aholder,iAii;
28: MatFactorInfo info;
29: PetscScalar *xsurf,*xint;
30: const PetscScalar *rxint;
31: #if defined(PETSC_USE_DEBUG_foo)
32: PetscScalar tmp;
33: #endif
34: PetscTable ht;
36: DMDAGetInfo(da,&dim,NULL,NULL,NULL,&mp,&np,&pp,&dof,NULL,NULL,NULL,NULL,NULL);
39: DMDAGetCorners(da,NULL,NULL,NULL,&m,&n,&p);
40: DMDAGetGhostCorners(da,&istart,&jstart,&kstart,&mwidth,&nwidth,&pwidth);
41: istart = istart ? -1 : 0;
42: jstart = jstart ? -1 : 0;
43: kstart = kstart ? -1 : 0;
45: /*
46: the columns of P are the interpolation of each coarse grid point (one for each vertex and edge)
47: to all the local degrees of freedom (this includes the vertices, edges and faces).
49: Xint are the subset of the interpolation into the interior
51: Xface are the interpolation onto faces but not into the interior
53: Xsurf are the interpolation onto the vertices and edges (the surfbasket)
54: Xint
55: Symbolically one could write P = (Xface) after interchanging the rows to match the natural ordering on the domain
56: Xsurf
57: */
58: N = (m - istart)*(n - jstart)*(p - kstart);
59: Nint = (m-2-istart)*(n-2-jstart)*(p-2-kstart);
60: Nface = 2*((m-2-istart)*(n-2-jstart) + (m-2-istart)*(p-2-kstart) + (n-2-jstart)*(p-2-kstart));
61: Nwire = 4*((m-2-istart) + (n-2-jstart) + (p-2-kstart)) + 8;
62: Nsurf = Nface + Nwire;
63: MatCreateSeqDense(MPI_COMM_SELF,Nint,26,NULL,&Xint);
64: MatCreateSeqDense(MPI_COMM_SELF,Nsurf,26,NULL,&Xsurf);
65: MatDenseGetArray(Xsurf,&xsurf);
67: /*
68: Require that all 12 edges and 6 faces have at least one grid point. Otherwise some of the columns of
69: Xsurf will be all zero (thus making the coarse matrix singular).
70: */
75: cnt = 0;
77: xsurf[cnt++] = 1;
78: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + Nsurf] = 1;
79: xsurf[cnt++ + 2*Nsurf] = 1;
81: for (j=1; j<n-1-jstart; j++) {
82: xsurf[cnt++ + 3*Nsurf] = 1;
83: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 4*Nsurf] = 1;
84: xsurf[cnt++ + 5*Nsurf] = 1;
85: }
87: xsurf[cnt++ + 6*Nsurf] = 1;
88: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 7*Nsurf] = 1;
89: xsurf[cnt++ + 8*Nsurf] = 1;
91: for (k=1; k<p-1-kstart; k++) {
92: xsurf[cnt++ + 9*Nsurf] = 1;
93: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 10*Nsurf] = 1;
94: xsurf[cnt++ + 11*Nsurf] = 1;
96: for (j=1; j<n-1-jstart; j++) {
97: xsurf[cnt++ + 12*Nsurf] = 1;
98: /* these are the interior nodes */
99: xsurf[cnt++ + 13*Nsurf] = 1;
100: }
102: xsurf[cnt++ + 14*Nsurf] = 1;
103: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 15*Nsurf] = 1;
104: xsurf[cnt++ + 16*Nsurf] = 1;
105: }
107: xsurf[cnt++ + 17*Nsurf] = 1;
108: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 18*Nsurf] = 1;
109: xsurf[cnt++ + 19*Nsurf] = 1;
111: for (j=1;j<n-1-jstart;j++) {
112: xsurf[cnt++ + 20*Nsurf] = 1;
113: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 21*Nsurf] = 1;
114: xsurf[cnt++ + 22*Nsurf] = 1;
115: }
117: xsurf[cnt++ + 23*Nsurf] = 1;
118: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 24*Nsurf] = 1;
119: xsurf[cnt++ + 25*Nsurf] = 1;
121: /* interpolations only sum to 1 when using direct solver */
122: #if defined(PETSC_USE_DEBUG_foo)
123: for (i=0; i<Nsurf; i++) {
124: tmp = 0.0;
125: for (j=0; j<26; j++) tmp += xsurf[i+j*Nsurf];
127: }
128: #endif
129: MatDenseRestoreArray(Xsurf,&xsurf);
130: /* MatView(Xsurf,PETSC_VIEWER_STDOUT_WORLD);*/
132: /*
133: I are the indices for all the needed vertices (in global numbering)
134: Iint are the indices for the interior values, I surf for the surface values
135: (This is just for the part of the global matrix obtained with MatCreateSubMatrix(), it
136: is NOT the local DMDA ordering.)
137: IIint and IIsurf are the same as the Iint, Isurf except they are in the global numbering
138: */
139: #define Endpoint(a,start,b) (a == 0 || a == (b-1-start))
140: PetscMalloc3(N,&II,Nint,&Iint,Nsurf,&Isurf);
141: PetscMalloc2(Nint,&IIint,Nsurf,&IIsurf);
142: for (k=0; k<p-kstart; k++) {
143: for (j=0; j<n-jstart; j++) {
144: for (i=0; i<m-istart; i++) {
145: II[c++] = i + j*mwidth + k*mwidth*nwidth;
147: if (!Endpoint(i,istart,m) && !Endpoint(j,jstart,n) && !Endpoint(k,kstart,p)) {
148: IIint[cint] = i + j*mwidth + k*mwidth*nwidth;
149: Iint[cint++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
150: } else {
151: IIsurf[csurf] = i + j*mwidth + k*mwidth*nwidth;
152: Isurf[csurf++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
153: }
154: }
155: }
156: }
160: DMGetLocalToGlobalMapping(da,<g);
161: ISLocalToGlobalMappingApply(ltg,N,II,II);
162: ISLocalToGlobalMappingApply(ltg,Nint,IIint,IIint);
163: ISLocalToGlobalMappingApply(ltg,Nsurf,IIsurf,IIsurf);
164: PetscObjectGetComm((PetscObject)da,&comm);
165: ISCreateGeneral(comm,N,II,PETSC_COPY_VALUES,&is);
166: ISCreateGeneral(PETSC_COMM_SELF,Nint,Iint,PETSC_COPY_VALUES,&isint);
167: ISCreateGeneral(PETSC_COMM_SELF,Nsurf,Isurf,PETSC_COPY_VALUES,&issurf);
168: PetscFree3(II,Iint,Isurf);
170: MatCreateSubMatrices(Aglobal,1,&is,&is,MAT_INITIAL_MATRIX,&Aholder);
171: A = *Aholder;
172: PetscFree(Aholder);
174: MatCreateSubMatrix(A,isint,isint,MAT_INITIAL_MATRIX,&Aii);
175: MatCreateSubMatrix(A,isint,issurf,MAT_INITIAL_MATRIX,&Ais);
176: MatCreateSubMatrix(A,issurf,isint,MAT_INITIAL_MATRIX,&Asi);
178: /*
179: Solve for the interpolation onto the interior Xint
180: */
181: MatMatMult(Ais,Xsurf,MAT_INITIAL_MATRIX,PETSC_DETERMINE,&Xint_tmp);
182: MatScale(Xint_tmp,-1.0);
183: if (exotic->directSolve) {
184: MatGetFactor(Aii,MATSOLVERPETSC,MAT_FACTOR_LU,&iAii);
185: MatFactorInfoInitialize(&info);
186: MatGetOrdering(Aii,MATORDERINGND,&row,&col);
187: MatLUFactorSymbolic(iAii,Aii,row,col,&info);
188: ISDestroy(&row);
189: ISDestroy(&col);
190: MatLUFactorNumeric(iAii,Aii,&info);
191: MatMatSolve(iAii,Xint_tmp,Xint);
192: MatDestroy(&iAii);
193: } else {
194: Vec b,x;
195: PetscScalar *xint_tmp;
197: MatDenseGetArray(Xint,&xint);
198: VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,NULL,&x);
199: MatDenseGetArray(Xint_tmp,&xint_tmp);
200: VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,NULL,&b);
201: KSPSetOperators(exotic->ksp,Aii,Aii);
202: for (i=0; i<26; i++) {
203: VecPlaceArray(x,xint+i*Nint);
204: VecPlaceArray(b,xint_tmp+i*Nint);
205: KSPSolve(exotic->ksp,b,x);
206: KSPCheckSolve(exotic->ksp,pc,x);
207: VecResetArray(x);
208: VecResetArray(b);
209: }
210: MatDenseRestoreArray(Xint,&xint);
211: MatDenseRestoreArray(Xint_tmp,&xint_tmp);
212: VecDestroy(&x);
213: VecDestroy(&b);
214: }
215: MatDestroy(&Xint_tmp);
217: #if defined(PETSC_USE_DEBUG_foo)
218: MatDenseGetArrayRead(Xint,&rxint);
219: for (i=0; i<Nint; i++) {
220: tmp = 0.0;
221: for (j=0; j<26; j++) tmp += rxint[i+j*Nint];
224: }
225: MatDenseRestoreArrayRead(Xint,&rxint);
226: /* MatView(Xint,PETSC_VIEWER_STDOUT_WORLD); */
227: #endif
229: /* total vertices total faces total edges */
230: Ntotal = (mp + 1)*(np + 1)*(pp + 1) + mp*np*(pp+1) + mp*pp*(np+1) + np*pp*(mp+1) + mp*(np+1)*(pp+1) + np*(mp+1)*(pp+1) + pp*(mp+1)*(np+1);
232: /*
233: For each vertex, edge, face on process (in the same orderings as used above) determine its local number including ghost points
234: */
235: cnt = 0;
237: gl[cnt++] = 0; { gl[cnt++] = 1;} gl[cnt++] = m-istart-1;
238: { gl[cnt++] = mwidth; { gl[cnt++] = mwidth+1;} gl[cnt++] = mwidth + m-istart-1;}
239: gl[cnt++] = mwidth*(n-jstart-1); { gl[cnt++] = mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*(n-jstart-1) + m-istart-1;
240: {
241: gl[cnt++] = mwidth*nwidth; { gl[cnt++] = mwidth*nwidth+1;} gl[cnt++] = mwidth*nwidth+ m-istart-1;
242: { gl[cnt++] = mwidth*nwidth + mwidth; /* these are the interior nodes */ gl[cnt++] = mwidth*nwidth + mwidth+m-istart-1;}
243: gl[cnt++] = mwidth*nwidth+ mwidth*(n-jstart-1); { gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1) + m-istart-1;
244: }
245: gl[cnt++] = mwidth*nwidth*(p-kstart-1); { gl[cnt++] = mwidth*nwidth*(p-kstart-1)+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1) + m-istart-1;
246: { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth; { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1)+mwidth+m-istart-1;}
247: gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth*(n-jstart-1); { gl[cnt++] = mwidth*nwidth*(p-kstart-1)+ mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth*(n-jstart-1) + m-istart-1;
249: /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */
250: /* convert that to global numbering and get them on all processes */
251: ISLocalToGlobalMappingApply(ltg,26,gl,gl);
252: /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */
253: PetscMalloc1(26*mp*np*pp,&globals);
254: MPI_Allgather(gl,26,MPIU_INT,globals,26,MPIU_INT,PetscObjectComm((PetscObject)da));
256: /* Number the coarse grid points from 0 to Ntotal */
257: MatGetSize(Aglobal,&Nt,NULL);
258: PetscTableCreate(Ntotal/3,Nt+1,&ht);
259: for (i=0; i<26*mp*np*pp; i++) {
260: PetscTableAddCount(ht,globals[i]+1);
261: }
262: PetscTableGetCount(ht,&cnt);
264: PetscFree(globals);
265: for (i=0; i<26; i++) {
266: PetscTableFind(ht,gl[i]+1,&gl[i]);
267: gl[i]--;
268: }
269: PetscTableDestroy(&ht);
270: /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */
272: /* construct global interpolation matrix */
273: MatGetLocalSize(Aglobal,&Ng,NULL);
274: if (reuse == MAT_INITIAL_MATRIX) {
275: MatCreateAIJ(PetscObjectComm((PetscObject)da),Ng,PETSC_DECIDE,PETSC_DECIDE,Ntotal,Nint+Nsurf,NULL,Nint+Nsurf,NULL,P);
276: } else {
277: MatZeroEntries(*P);
278: }
279: MatSetOption(*P,MAT_ROW_ORIENTED,PETSC_FALSE);
280: MatDenseGetArrayRead(Xint,&rxint);
281: MatSetValues(*P,Nint,IIint,26,gl,rxint,INSERT_VALUES);
282: MatDenseRestoreArrayRead(Xint,&rxint);
283: MatDenseGetArrayRead(Xsurf,&rxint);
284: MatSetValues(*P,Nsurf,IIsurf,26,gl,rxint,INSERT_VALUES);
285: MatDenseRestoreArrayRead(Xsurf,&rxint);
286: MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);
287: MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);
288: PetscFree2(IIint,IIsurf);
290: #if defined(PETSC_USE_DEBUG_foo)
291: {
292: Vec x,y;
293: PetscScalar *yy;
294: VecCreateMPI(PetscObjectComm((PetscObject)da),Ng,PETSC_DETERMINE,&y);
295: VecCreateMPI(PetscObjectComm((PetscObject)da),PETSC_DETERMINE,Ntotal,&x);
296: VecSet(x,1.0);
297: MatMult(*P,x,y);
298: VecGetArray(y,&yy);
299: for (i=0; i<Ng; i++) {
301: }
302: VecRestoreArray(y,&yy);
303: VecDestroy(x);
304: VecDestroy(y);
305: }
306: #endif
308: MatDestroy(&Aii);
309: MatDestroy(&Ais);
310: MatDestroy(&Asi);
311: MatDestroy(&A);
312: ISDestroy(&is);
313: ISDestroy(&isint);
314: ISDestroy(&issurf);
315: MatDestroy(&Xint);
316: MatDestroy(&Xsurf);
317: return 0;
318: }
320: /*
321: DMDAGetFaceInterpolation - Gets the interpolation for a face based coarse space
323: */
324: PetscErrorCode DMDAGetFaceInterpolation(PC pc,DM da,PC_Exotic *exotic,Mat Aglobal,MatReuse reuse,Mat *P)
325: {
326: PetscInt dim,i,j,k,m,n,p,dof,Nint,Nface,Nwire,Nsurf,*Iint,*Isurf,cint = 0,csurf = 0,istart,jstart,kstart,*II,N,c = 0;
327: PetscInt mwidth,nwidth,pwidth,cnt,mp,np,pp,Ntotal,gl[6],*globals,Ng,*IIint,*IIsurf,Nt;
328: Mat Xint, Xsurf,Xint_tmp;
329: IS isint,issurf,is,row,col;
330: ISLocalToGlobalMapping ltg;
331: MPI_Comm comm;
332: Mat A,Aii,Ais,Asi,*Aholder,iAii;
333: MatFactorInfo info;
334: PetscScalar *xsurf,*xint;
335: const PetscScalar *rxint;
336: #if defined(PETSC_USE_DEBUG_foo)
337: PetscScalar tmp;
338: #endif
339: PetscTable ht;
341: DMDAGetInfo(da,&dim,NULL,NULL,NULL,&mp,&np,&pp,&dof,NULL,NULL,NULL,NULL,NULL);
344: DMDAGetCorners(da,NULL,NULL,NULL,&m,&n,&p);
345: DMDAGetGhostCorners(da,&istart,&jstart,&kstart,&mwidth,&nwidth,&pwidth);
346: istart = istart ? -1 : 0;
347: jstart = jstart ? -1 : 0;
348: kstart = kstart ? -1 : 0;
350: /*
351: the columns of P are the interpolation of each coarse grid point (one for each vertex and edge)
352: to all the local degrees of freedom (this includes the vertices, edges and faces).
354: Xint are the subset of the interpolation into the interior
356: Xface are the interpolation onto faces but not into the interior
358: Xsurf are the interpolation onto the vertices and edges (the surfbasket)
359: Xint
360: Symbolically one could write P = (Xface) after interchanging the rows to match the natural ordering on the domain
361: Xsurf
362: */
363: N = (m - istart)*(n - jstart)*(p - kstart);
364: Nint = (m-2-istart)*(n-2-jstart)*(p-2-kstart);
365: Nface = 2*((m-2-istart)*(n-2-jstart) + (m-2-istart)*(p-2-kstart) + (n-2-jstart)*(p-2-kstart));
366: Nwire = 4*((m-2-istart) + (n-2-jstart) + (p-2-kstart)) + 8;
367: Nsurf = Nface + Nwire;
368: MatCreateSeqDense(MPI_COMM_SELF,Nint,6,NULL,&Xint);
369: MatCreateSeqDense(MPI_COMM_SELF,Nsurf,6,NULL,&Xsurf);
370: MatDenseGetArray(Xsurf,&xsurf);
372: /*
373: Require that all 12 edges and 6 faces have at least one grid point. Otherwise some of the columns of
374: Xsurf will be all zero (thus making the coarse matrix singular).
375: */
380: cnt = 0;
381: for (j=1; j<n-1-jstart; j++) {
382: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 0*Nsurf] = 1;
383: }
385: for (k=1; k<p-1-kstart; k++) {
386: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 1*Nsurf] = 1;
387: for (j=1; j<n-1-jstart; j++) {
388: xsurf[cnt++ + 2*Nsurf] = 1;
389: /* these are the interior nodes */
390: xsurf[cnt++ + 3*Nsurf] = 1;
391: }
392: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 4*Nsurf] = 1;
393: }
394: for (j=1;j<n-1-jstart;j++) {
395: for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 5*Nsurf] = 1;
396: }
398: #if defined(PETSC_USE_DEBUG_foo)
399: for (i=0; i<Nsurf; i++) {
400: tmp = 0.0;
401: for (j=0; j<6; j++) tmp += xsurf[i+j*Nsurf];
404: }
405: #endif
406: MatDenseRestoreArray(Xsurf,&xsurf);
407: /* MatView(Xsurf,PETSC_VIEWER_STDOUT_WORLD);*/
409: /*
410: I are the indices for all the needed vertices (in global numbering)
411: Iint are the indices for the interior values, I surf for the surface values
412: (This is just for the part of the global matrix obtained with MatCreateSubMatrix(), it
413: is NOT the local DMDA ordering.)
414: IIint and IIsurf are the same as the Iint, Isurf except they are in the global numbering
415: */
416: #define Endpoint(a,start,b) (a == 0 || a == (b-1-start))
417: PetscMalloc3(N,&II,Nint,&Iint,Nsurf,&Isurf);
418: PetscMalloc2(Nint,&IIint,Nsurf,&IIsurf);
419: for (k=0; k<p-kstart; k++) {
420: for (j=0; j<n-jstart; j++) {
421: for (i=0; i<m-istart; i++) {
422: II[c++] = i + j*mwidth + k*mwidth*nwidth;
424: if (!Endpoint(i,istart,m) && !Endpoint(j,jstart,n) && !Endpoint(k,kstart,p)) {
425: IIint[cint] = i + j*mwidth + k*mwidth*nwidth;
426: Iint[cint++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
427: } else {
428: IIsurf[csurf] = i + j*mwidth + k*mwidth*nwidth;
429: Isurf[csurf++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
430: }
431: }
432: }
433: }
437: DMGetLocalToGlobalMapping(da,<g);
438: ISLocalToGlobalMappingApply(ltg,N,II,II);
439: ISLocalToGlobalMappingApply(ltg,Nint,IIint,IIint);
440: ISLocalToGlobalMappingApply(ltg,Nsurf,IIsurf,IIsurf);
441: PetscObjectGetComm((PetscObject)da,&comm);
442: ISCreateGeneral(comm,N,II,PETSC_COPY_VALUES,&is);
443: ISCreateGeneral(PETSC_COMM_SELF,Nint,Iint,PETSC_COPY_VALUES,&isint);
444: ISCreateGeneral(PETSC_COMM_SELF,Nsurf,Isurf,PETSC_COPY_VALUES,&issurf);
445: PetscFree3(II,Iint,Isurf);
447: ISSort(is);
448: MatCreateSubMatrices(Aglobal,1,&is,&is,MAT_INITIAL_MATRIX,&Aholder);
449: A = *Aholder;
450: PetscFree(Aholder);
452: MatCreateSubMatrix(A,isint,isint,MAT_INITIAL_MATRIX,&Aii);
453: MatCreateSubMatrix(A,isint,issurf,MAT_INITIAL_MATRIX,&Ais);
454: MatCreateSubMatrix(A,issurf,isint,MAT_INITIAL_MATRIX,&Asi);
456: /*
457: Solve for the interpolation onto the interior Xint
458: */
459: MatMatMult(Ais,Xsurf,MAT_INITIAL_MATRIX,PETSC_DETERMINE,&Xint_tmp);
460: MatScale(Xint_tmp,-1.0);
462: if (exotic->directSolve) {
463: MatGetFactor(Aii,MATSOLVERPETSC,MAT_FACTOR_LU,&iAii);
464: MatFactorInfoInitialize(&info);
465: MatGetOrdering(Aii,MATORDERINGND,&row,&col);
466: MatLUFactorSymbolic(iAii,Aii,row,col,&info);
467: ISDestroy(&row);
468: ISDestroy(&col);
469: MatLUFactorNumeric(iAii,Aii,&info);
470: MatMatSolve(iAii,Xint_tmp,Xint);
471: MatDestroy(&iAii);
472: } else {
473: Vec b,x;
474: PetscScalar *xint_tmp;
476: MatDenseGetArray(Xint,&xint);
477: VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,NULL,&x);
478: MatDenseGetArray(Xint_tmp,&xint_tmp);
479: VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,NULL,&b);
480: KSPSetOperators(exotic->ksp,Aii,Aii);
481: for (i=0; i<6; i++) {
482: VecPlaceArray(x,xint+i*Nint);
483: VecPlaceArray(b,xint_tmp+i*Nint);
484: KSPSolve(exotic->ksp,b,x);
485: KSPCheckSolve(exotic->ksp,pc,x);
486: VecResetArray(x);
487: VecResetArray(b);
488: }
489: MatDenseRestoreArray(Xint,&xint);
490: MatDenseRestoreArray(Xint_tmp,&xint_tmp);
491: VecDestroy(&x);
492: VecDestroy(&b);
493: }
494: MatDestroy(&Xint_tmp);
496: #if defined(PETSC_USE_DEBUG_foo)
497: MatDenseGetArrayRead(Xint,&rxint);
498: for (i=0; i<Nint; i++) {
499: tmp = 0.0;
500: for (j=0; j<6; j++) tmp += rxint[i+j*Nint];
503: }
504: MatDenseRestoreArrayRead(Xint,&rxint);
505: /* MatView(Xint,PETSC_VIEWER_STDOUT_WORLD); */
506: #endif
508: /* total faces */
509: Ntotal = mp*np*(pp+1) + mp*pp*(np+1) + np*pp*(mp+1);
511: /*
512: For each vertex, edge, face on process (in the same orderings as used above) determine its local number including ghost points
513: */
514: cnt = 0;
515: { gl[cnt++] = mwidth+1;}
516: {
517: { gl[cnt++] = mwidth*nwidth+1;}
518: { gl[cnt++] = mwidth*nwidth + mwidth; /* these are the interior nodes */ gl[cnt++] = mwidth*nwidth + mwidth+m-istart-1;}
519: { gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1)+1;}
520: }
521: { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth+1;}
523: /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */
524: /* convert that to global numbering and get them on all processes */
525: ISLocalToGlobalMappingApply(ltg,6,gl,gl);
526: /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */
527: PetscMalloc1(6*mp*np*pp,&globals);
528: MPI_Allgather(gl,6,MPIU_INT,globals,6,MPIU_INT,PetscObjectComm((PetscObject)da));
530: /* Number the coarse grid points from 0 to Ntotal */
531: MatGetSize(Aglobal,&Nt,NULL);
532: PetscTableCreate(Ntotal/3,Nt+1,&ht);
533: for (i=0; i<6*mp*np*pp; i++) {
534: PetscTableAddCount(ht,globals[i]+1);
535: }
536: PetscTableGetCount(ht,&cnt);
538: PetscFree(globals);
539: for (i=0; i<6; i++) {
540: PetscTableFind(ht,gl[i]+1,&gl[i]);
541: gl[i]--;
542: }
543: PetscTableDestroy(&ht);
544: /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */
546: /* construct global interpolation matrix */
547: MatGetLocalSize(Aglobal,&Ng,NULL);
548: if (reuse == MAT_INITIAL_MATRIX) {
549: MatCreateAIJ(PetscObjectComm((PetscObject)da),Ng,PETSC_DECIDE,PETSC_DECIDE,Ntotal,Nint+Nsurf,NULL,Nint,NULL,P);
550: } else {
551: MatZeroEntries(*P);
552: }
553: MatSetOption(*P,MAT_ROW_ORIENTED,PETSC_FALSE);
554: MatDenseGetArrayRead(Xint,&rxint);
555: MatSetValues(*P,Nint,IIint,6,gl,rxint,INSERT_VALUES);
556: MatDenseRestoreArrayRead(Xint,&rxint);
557: MatDenseGetArrayRead(Xsurf,&rxint);
558: MatSetValues(*P,Nsurf,IIsurf,6,gl,rxint,INSERT_VALUES);
559: MatDenseRestoreArrayRead(Xsurf,&rxint);
560: MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);
561: MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);
562: PetscFree2(IIint,IIsurf);
564: #if defined(PETSC_USE_DEBUG_foo)
565: {
566: Vec x,y;
567: PetscScalar *yy;
568: VecCreateMPI(PetscObjectComm((PetscObject)da),Ng,PETSC_DETERMINE,&y);
569: VecCreateMPI(PetscObjectComm((PetscObject)da),PETSC_DETERMINE,Ntotal,&x);
570: VecSet(x,1.0);
571: MatMult(*P,x,y);
572: VecGetArray(y,&yy);
573: for (i=0; i<Ng; i++) {
575: }
576: VecRestoreArray(y,&yy);
577: VecDestroy(x);
578: VecDestroy(y);
579: }
580: #endif
582: MatDestroy(&Aii);
583: MatDestroy(&Ais);
584: MatDestroy(&Asi);
585: MatDestroy(&A);
586: ISDestroy(&is);
587: ISDestroy(&isint);
588: ISDestroy(&issurf);
589: MatDestroy(&Xint);
590: MatDestroy(&Xsurf);
591: return 0;
592: }
594: /*@
595: PCExoticSetType - Sets the type of coarse grid interpolation to use
597: Logically Collective on PC
599: Input Parameters:
600: + pc - the preconditioner context
601: - type - either PC_EXOTIC_FACE or PC_EXOTIC_WIREBASKET (defaults to face)
603: Notes:
604: The face based interpolation has 1 degree of freedom per face and ignores the
605: edge and vertex values completely in the coarse problem. For any seven point
606: stencil the interpolation of a constant on all faces into the interior is that constant.
608: The wirebasket interpolation has 1 degree of freedom per vertex, per edge and
609: per face. A constant on the subdomain boundary is interpolated as that constant
610: in the interior of the domain.
612: The coarse grid matrix is obtained via the Galerkin computation A_c = R A R^T, hence
613: if A is nonsingular A_c is also nonsingular.
615: Both interpolations are suitable for only scalar problems.
617: Level: intermediate
619: .seealso: PCEXOTIC, PCExoticType()
620: @*/
621: PetscErrorCode PCExoticSetType(PC pc,PCExoticType type)
622: {
625: PetscTryMethod(pc,"PCExoticSetType_C",(PC,PCExoticType),(pc,type));
626: return 0;
627: }
629: static PetscErrorCode PCExoticSetType_Exotic(PC pc,PCExoticType type)
630: {
631: PC_MG *mg = (PC_MG*)pc->data;
632: PC_Exotic *ctx = (PC_Exotic*) mg->innerctx;
634: ctx->type = type;
635: return 0;
636: }
638: PetscErrorCode PCSetUp_Exotic(PC pc)
639: {
640: Mat A;
641: PC_MG *mg = (PC_MG*)pc->data;
642: PC_Exotic *ex = (PC_Exotic*) mg->innerctx;
643: MatReuse reuse = (ex->P) ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
646: PCGetOperators(pc,NULL,&A);
647: if (ex->type == PC_EXOTIC_FACE) {
648: DMDAGetFaceInterpolation(pc,pc->dm,ex,A,reuse,&ex->P);
649: } else if (ex->type == PC_EXOTIC_WIREBASKET) {
650: DMDAGetWireBasketInterpolation(pc,pc->dm,ex,A,reuse,&ex->P);
651: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unknown exotic coarse space %d",ex->type);
652: PCMGSetInterpolation(pc,1,ex->P);
653: /* if PC has attached DM we must remove it or the PCMG will use it to compute incorrect sized vectors and interpolations */
654: PCSetDM(pc,NULL);
655: PCSetUp_MG(pc);
656: return 0;
657: }
659: PetscErrorCode PCDestroy_Exotic(PC pc)
660: {
661: PC_MG *mg = (PC_MG*)pc->data;
662: PC_Exotic *ctx = (PC_Exotic*) mg->innerctx;
664: MatDestroy(&ctx->P);
665: KSPDestroy(&ctx->ksp);
666: PetscFree(ctx);
667: PCDestroy_MG(pc);
668: return 0;
669: }
671: PetscErrorCode PCView_Exotic(PC pc,PetscViewer viewer)
672: {
673: PC_MG *mg = (PC_MG*)pc->data;
674: PetscBool iascii;
675: PC_Exotic *ctx = (PC_Exotic*) mg->innerctx;
677: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
678: if (iascii) {
679: PetscViewerASCIIPrintf(viewer," Exotic type = %s\n",PCExoticTypes[ctx->type]);
680: if (ctx->directSolve) {
681: PetscViewerASCIIPrintf(viewer," Using direct solver to construct interpolation\n");
682: } else {
683: PetscViewer sviewer;
684: PetscMPIInt rank;
686: PetscViewerASCIIPrintf(viewer," Using iterative solver to construct interpolation\n");
687: PetscViewerASCIIPushTab(viewer);
688: PetscViewerASCIIPushTab(viewer); /* should not need to push this twice? */
689: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
690: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
691: if (rank == 0) {
692: KSPView(ctx->ksp,sviewer);
693: }
694: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
695: PetscViewerASCIIPopTab(viewer);
696: PetscViewerASCIIPopTab(viewer);
697: }
698: }
699: PCView_MG(pc,viewer);
700: return 0;
701: }
703: PetscErrorCode PCSetFromOptions_Exotic(PetscOptionItems *PetscOptionsObject,PC pc)
704: {
705: PetscBool flg;
706: PC_MG *mg = (PC_MG*)pc->data;
707: PCExoticType mgctype;
708: PC_Exotic *ctx = (PC_Exotic*) mg->innerctx;
710: PetscOptionsHead(PetscOptionsObject,"Exotic coarse space options");
711: PetscOptionsEnum("-pc_exotic_type","face or wirebasket","PCExoticSetType",PCExoticTypes,(PetscEnum)ctx->type,(PetscEnum*)&mgctype,&flg);
712: if (flg) {
713: PCExoticSetType(pc,mgctype);
714: }
715: PetscOptionsBool("-pc_exotic_direct_solver","use direct solver to construct interpolation","None",ctx->directSolve,&ctx->directSolve,NULL);
716: if (!ctx->directSolve) {
717: if (!ctx->ksp) {
718: const char *prefix;
719: KSPCreate(PETSC_COMM_SELF,&ctx->ksp);
720: KSPSetErrorIfNotConverged(ctx->ksp,pc->erroriffailure);
721: PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)pc,1);
722: PetscLogObjectParent((PetscObject)pc,(PetscObject)ctx->ksp);
723: PCGetOptionsPrefix(pc,&prefix);
724: KSPSetOptionsPrefix(ctx->ksp,prefix);
725: KSPAppendOptionsPrefix(ctx->ksp,"exotic_");
726: }
727: KSPSetFromOptions(ctx->ksp);
728: }
729: PetscOptionsTail();
730: return 0;
731: }
733: /*MC
734: PCEXOTIC - Two level overlapping Schwarz preconditioner with exotic (non-standard) coarse grid spaces
736: This uses the PCMG infrastructure restricted to two levels and the face and wirebasket based coarse
737: grid spaces.
739: Notes:
740: By default this uses GMRES on the fine grid smoother so this should be used with KSPFGMRES or the smoother changed to not use GMRES
742: References:
743: + * - These coarse grid spaces originate in the work of Bramble, Pasciak and Schatz, "The Construction
744: of Preconditioners for Elliptic Problems by Substructing IV", Mathematics of Computation, volume 53, 1989.
745: . * - They were generalized slightly in "Domain Decomposition Method for Linear Elasticity", Ph. D. thesis, Barry Smith,
746: New York University, 1990.
747: . * - They were then explored in great detail in Dryja, Smith, Widlund, "Schwarz Analysis
748: of Iterative Substructuring Methods for Elliptic Problems in Three Dimensions, SIAM Journal on Numerical
749: Analysis, volume 31. 1994. These were developed in the context of iterative substructuring preconditioners.
750: . * - They were then ingeniously applied as coarse grid spaces for overlapping Schwarz methods by Dohrmann and Widlund.
751: They refer to them as GDSW (generalized Dryja, Smith, Widlund preconditioners). See, for example,
752: Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. Extending theory for domain decomposition algorithms to irregular subdomains. In Ulrich Langer, Marco
753: Discacciati, David Keyes, Olof Widlund, and Walter Zulehner, editors, Proceedings
754: of the 17th International Conference on Domain Decomposition Methods in
755: Science and Engineering, held in Strobl, Austria, 2006, number 60 in
756: Springer Verlag, Lecture Notes in Computational Science and Engineering, 2007.
757: . * - Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. A family of energy minimizing coarse spaces for overlapping Schwarz preconditioners. In Ulrich Langer,
758: Marco Discacciati, David Keyes, Olof Widlund, and Walter Zulehner, editors, Proceedings
759: of the 17th International Conference on Domain Decomposition Methods
760: in Science and Engineering, held in Strobl, Austria, 2006, number 60 in
761: Springer Verlag, Lecture Notes in Computational Science and Engineering, 2007
762: . * - Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. Domain decomposition
763: for less regular subdomains: Overlapping Schwarz in two dimensions. SIAM J.
764: Numer. Anal., 46(4), 2008.
765: - * - Clark R. Dohrmann and Olof B. Widlund. An overlapping Schwarz
766: algorithm for almost incompressible elasticity. Technical Report
767: TR2008 912, Department of Computer Science, Courant Institute
768: of Mathematical Sciences, New York University, May 2008. URL:
770: Options Database: The usual PCMG options are supported, such as -mg_levels_pc_type <type> -mg_coarse_pc_type <type>
771: -pc_mg_type <type>
773: Level: advanced
775: .seealso: PCMG, PCSetDM(), PCExoticType, PCExoticSetType()
776: M*/
778: PETSC_EXTERN PetscErrorCode PCCreate_Exotic(PC pc)
779: {
780: PC_Exotic *ex;
781: PC_MG *mg;
783: /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */
784: if (pc->ops->destroy) {
785: (*pc->ops->destroy)(pc);
786: pc->data = NULL;
787: }
788: PetscFree(((PetscObject)pc)->type_name);
789: ((PetscObject)pc)->type_name = NULL;
791: PCSetType(pc,PCMG);
792: PCMGSetLevels(pc,2,NULL);
793: PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
794: PetscNew(&ex); \
795: ex->type = PC_EXOTIC_FACE;
796: mg = (PC_MG*) pc->data;
797: mg->innerctx = ex;
799: pc->ops->setfromoptions = PCSetFromOptions_Exotic;
800: pc->ops->view = PCView_Exotic;
801: pc->ops->destroy = PCDestroy_Exotic;
802: pc->ops->setup = PCSetUp_Exotic;
804: PetscObjectComposeFunction((PetscObject)pc,"PCExoticSetType_C",PCExoticSetType_Exotic);
805: return 0;
806: }