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