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