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