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, &ltg);
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, &ltg);
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: }