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