Actual source code: ex73.c


  2: /*
  3: This example was derived from src/ksp/ksp/tutorials ex29.c

  5: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation

  7:    -div \rho grad u = f,  0 < x,y < 1,

  9: with forcing function

 11:    f = e^{-x^2/\nu} e^{-y^2/\nu}

 13: with Dirichlet boundary conditions

 15:    u = f(x,y) for x = 0, x = 1, y = 0, y = 1

 17: or pure Neumman boundary conditions.
 18: */

 20: static char help[] = "Solves 2D inhomogeneous Laplacian. Demonstates using PCTelescopeSetCoarseDM functionality of PCTelescope via a DMShell\n\n";

 22: #include <petscdm.h>
 23: #include <petscdmda.h>
 24: #include <petscdmshell.h>
 25: #include <petscksp.h>

 27: PetscErrorCode ComputeMatrix_ShellDA(KSP,Mat,Mat,void*);
 28: PetscErrorCode ComputeMatrix_DMDA(DM,Mat,Mat,void*);
 29: PetscErrorCode ComputeRHS_DMDA(DM,Vec,void*);
 30: PetscErrorCode DMShellCreate_ShellDA(DM,DM*);
 31: PetscErrorCode DMFieldScatter_ShellDA(DM,Vec,ScatterMode,DM,Vec);
 32: PetscErrorCode DMStateScatter_ShellDA(DM,ScatterMode,DM);

 34: typedef enum {DIRICHLET, NEUMANN} BCType;

 36: typedef struct {
 37:   PetscReal rho;
 38:   PetscReal nu;
 39:   BCType    bcType;
 40:   MPI_Comm  comm;
 41: } UserContext;

 43: PetscErrorCode UserContextCreate(MPI_Comm comm,UserContext **ctx)
 44: {
 45:   UserContext    *user;
 46:   const char     *bcTypes[2] = {"dirichlet","neumann"};
 47:   PetscInt       bc;

 51:   PetscCalloc1(1,&user);
 52:   user->comm = comm;
 53:   PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 54:   user->rho = 1.0;
 55:   PetscOptionsReal("-rho", "The conductivity", "ex29.c", user->rho, &user->rho, NULL);
 56:   user->nu = 0.1;
 57:   PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user->nu, &user->nu, NULL);
 58:   bc = (PetscInt)DIRICHLET;
 59:   PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);
 60:   user->bcType = (BCType)bc;
 61:   PetscOptionsEnd();
 62:   *ctx = user;
 63:   return 0;
 64: }

 66: PetscErrorCode CommCoarsen(MPI_Comm comm,PetscInt number,PetscSubcomm *p)
 67: {
 68:   PetscSubcomm   psubcomm;
 70:   PetscSubcommCreate(comm,&psubcomm);
 71:   PetscSubcommSetNumber(psubcomm,number);
 72:   PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_INTERLACED);
 73:   *p = psubcomm;
 74:   return 0;
 75: }

 77: PetscErrorCode CommHierarchyCreate(MPI_Comm comm,PetscInt n,PetscInt number[],PetscSubcomm pscommlist[])
 78: {
 79:   PetscInt       k;
 80:   PetscBool      view_hierarchy = PETSC_FALSE;

 83:   for (k=0; k<n; k++) {
 84:     pscommlist[k] = NULL;
 85:   }

 87:   if (n < 1) return 0;

 89:   CommCoarsen(comm,number[n-1],&pscommlist[n-1]);
 90:   for (k=n-2; k>=0; k--) {
 91:     MPI_Comm comm_k = PetscSubcommChild(pscommlist[k+1]);
 92:     if (pscommlist[k+1]->color == 0) {
 93:       CommCoarsen(comm_k,number[k],&pscommlist[k]);
 94:     }
 95:   }

 97:   PetscOptionsGetBool(NULL,NULL,"-view_hierarchy",&view_hierarchy,NULL);
 98:   if (view_hierarchy) {
 99:     PetscMPIInt size;

101:     MPI_Comm_size(comm,&size);
102:     PetscPrintf(comm,"level[%D] size %d\n",n,(int)size);
103:     for (k=n-1; k>=0; k--) {
104:       if (pscommlist[k]) {
105:         MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);

107:         if (pscommlist[k]->color == 0) {
108:           MPI_Comm_size(comm_k,&size);
109:           PetscPrintf(comm_k,"level[%D] size %d\n",k,(int)size);
110:         }
111:       }
112:     }
113:   }
114:   return 0;
115: }

117: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
118: static PetscErrorCode _DMDADetermineRankFromGlobalIJ_2d(PetscInt i,PetscInt j,PetscInt Mp,PetscInt Np,
119:                                                         PetscInt start_i[],PetscInt start_j[],
120:                                                         PetscInt span_i[],PetscInt span_j[],
121:                                                         PetscMPIInt *_pi,PetscMPIInt *_pj,PetscMPIInt *rank_re)
122: {
123:   PetscInt pi,pj,n;

126:   *rank_re = -1;
127:   pi = pj = -1;
128:   if (_pi) {
129:     for (n=0; n<Mp; n++) {
130:       if ((i >= start_i[n]) && (i < start_i[n]+span_i[n])) {
131:         pi = n;
132:         break;
133:       }
134:     }
136:     *_pi = (PetscMPIInt)pi;
137:   }

139:   if (_pj) {
140:     for (n=0; n<Np; n++) {
141:       if ((j >= start_j[n]) && (j < start_j[n]+span_j[n])) {
142:         pj = n;
143:         break;
144:       }
145:     }
147:     *_pj = (PetscMPIInt)pj;
148:   }

150:   *rank_re = (PetscMPIInt)(pi + pj * Mp);
151:   return 0;
152: }

154: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
155: static PetscErrorCode _DMDADetermineGlobalS0_2d(PetscMPIInt rank_re,PetscInt Mp_re,PetscInt Np_re,
156:                                                 PetscInt range_i_re[],PetscInt range_j_re[],PetscInt *s0)
157: {
158:   PetscInt    i,j,start_IJ = 0;
159:   PetscMPIInt rank_ij;

162:   *s0 = -1;
163:   for (j=0; j<Np_re; j++) {
164:     for (i=0; i<Mp_re; i++) {
165:       rank_ij = (PetscMPIInt)(i + j*Mp_re);
166:       if (rank_ij < rank_re) {
167:         start_IJ += range_i_re[i]*range_j_re[j];
168:       }
169:     }
170:   }
171:   *s0 = start_IJ;
172:   return 0;
173: }

175: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
176: static PetscErrorCode DMDACreatePermutation_2d(DM dmrepart,DM dmf,Mat *mat)
177: {
179:   PetscInt       k,sum,Mp_re = 0,Np_re = 0;
180:   PetscInt       nx,ny,sr,er,Mr,ndof;
181:   PetscInt       i,j,location,startI[2],endI[2],lenI[2];
182:   const PetscInt *_range_i_re = NULL,*_range_j_re = NULL;
183:   PetscInt       *range_i_re,*range_j_re;
184:   PetscInt       *start_i_re,*start_j_re;
185:   MPI_Comm       comm;
186:   Vec            V;
187:   Mat            Pscalar;

190:   PetscInfo(dmf,"setting up the permutation matrix (DMDA-2D)\n");
191:   PetscObjectGetComm((PetscObject)dmf,&comm);

193:   _range_i_re = _range_j_re = NULL;
194:   /* Create DMDA on the child communicator */
195:   if (dmrepart) {
196:     DMDAGetInfo(dmrepart,NULL,NULL,NULL,NULL,&Mp_re,&Np_re,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
197:     DMDAGetOwnershipRanges(dmrepart,&_range_i_re,&_range_j_re,NULL);
198:   }

200:   /* note - assume rank 0 always participates */
201:   MPI_Bcast(&Mp_re,1,MPIU_INT,0,comm);
202:   MPI_Bcast(&Np_re,1,MPIU_INT,0,comm);

204:   PetscCalloc1(Mp_re,&range_i_re);
205:   PetscCalloc1(Np_re,&range_j_re);

207:   if (_range_i_re) PetscArraycpy(range_i_re,_range_i_re,Mp_re);
208:   if (_range_j_re) PetscArraycpy(range_j_re,_range_j_re,Np_re);

210:   MPI_Bcast(range_i_re,Mp_re,MPIU_INT,0,comm);
211:   MPI_Bcast(range_j_re,Np_re,MPIU_INT,0,comm);

213:   PetscMalloc1(Mp_re,&start_i_re);
214:   PetscMalloc1(Np_re,&start_j_re);

216:   sum = 0;
217:   for (k=0; k<Mp_re; k++) {
218:     start_i_re[k] = sum;
219:     sum += range_i_re[k];
220:   }

222:   sum = 0;
223:   for (k=0; k<Np_re; k++) {
224:     start_j_re[k] = sum;
225:     sum += range_j_re[k];
226:   }

228:   /* Create permutation */
229:   DMDAGetInfo(dmf,NULL,&nx,&ny,NULL,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);
230:   DMGetGlobalVector(dmf,&V);
231:   VecGetSize(V,&Mr);
232:   VecGetOwnershipRange(V,&sr,&er);
233:   DMRestoreGlobalVector(dmf,&V);
234:   sr = sr / ndof;
235:   er = er / ndof;
236:   Mr = Mr / ndof;

238:   MatCreate(comm,&Pscalar);
239:   MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);
240:   MatSetType(Pscalar,MATAIJ);
241:   MatSeqAIJSetPreallocation(Pscalar,1,NULL);
242:   MatMPIAIJSetPreallocation(Pscalar,1,NULL,1,NULL);

244:   DMDAGetCorners(dmf,NULL,NULL,NULL,&lenI[0],&lenI[1],NULL);
245:   DMDAGetCorners(dmf,&startI[0],&startI[1],NULL,&endI[0],&endI[1],NULL);
246:   endI[0] += startI[0];
247:   endI[1] += startI[1];

249:   for (j=startI[1]; j<endI[1]; j++) {
250:     for (i=startI[0]; i<endI[0]; i++) {
251:       PetscMPIInt rank_ijk_re,rank_reI[] = {0,0};
252:       PetscInt    s0_re;
253:       PetscInt    ii,jj,local_ijk_re,mapped_ijk;
254:       PetscInt    lenI_re[] = {0,0};

256:       location = (i - startI[0]) + (j - startI[1])*lenI[0];
257:       _DMDADetermineRankFromGlobalIJ_2d(i,j,Mp_re,Np_re,
258:                                              start_i_re,start_j_re,
259:                                              range_i_re,range_j_re,
260:                                              &rank_reI[0],&rank_reI[1],&rank_ijk_re);

262:       _DMDADetermineGlobalS0_2d(rank_ijk_re,Mp_re,Np_re,range_i_re,range_j_re,&s0_re);

264:       ii = i - start_i_re[ rank_reI[0] ];
266:       jj = j - start_j_re[ rank_reI[1] ];

269:       lenI_re[0] = range_i_re[ rank_reI[0] ];
270:       lenI_re[1] = range_j_re[ rank_reI[1] ];
271:       local_ijk_re = ii + jj * lenI_re[0];
272:       mapped_ijk = s0_re + local_ijk_re;
273:       MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);
274:     }
275:   }
276:   MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);
277:   MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);

279:   *mat = Pscalar;

281:   PetscFree(range_i_re);
282:   PetscFree(range_j_re);
283:   PetscFree(start_i_re);
284:   PetscFree(start_j_re);
285:   return 0;
286: }

288: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
289: static PetscErrorCode PCTelescopeSetUp_dmda_scatters(DM dmf,DM dmc)
290: {
291:   Vec            xred,yred,xtmp,x,xp;
292:   VecScatter     scatter;
293:   IS             isin;
294:   PetscInt       m,bs,st,ed;
295:   MPI_Comm       comm;
296:   VecType        vectype;

299:   PetscObjectGetComm((PetscObject)dmf,&comm);
300:   DMCreateGlobalVector(dmf,&x);
301:   VecGetBlockSize(x,&bs);
302:   VecGetType(x,&vectype);

304:   /* cannot use VecDuplicate as xp is already composed with dmf */
305:   /*VecDuplicate(x,&xp);*/
306:   {
307:     PetscInt m,M;

309:     VecGetSize(x,&M);
310:     VecGetLocalSize(x,&m);
311:     VecCreate(comm,&xp);
312:     VecSetSizes(xp,m,M);
313:     VecSetBlockSize(xp,bs);
314:     VecSetType(xp,vectype);
315:   }

317:   m = 0;
318:   xred = NULL;
319:   yred = NULL;
320:   if (dmc) {
321:     DMCreateGlobalVector(dmc,&xred);
322:     VecDuplicate(xred,&yred);
323:     VecGetOwnershipRange(xred,&st,&ed);
324:     ISCreateStride(comm,ed-st,st,1,&isin);
325:     VecGetLocalSize(xred,&m);
326:   } else {
327:     VecGetOwnershipRange(x,&st,&ed);
328:     ISCreateStride(comm,0,st,1,&isin);
329:   }
330:   ISSetBlockSize(isin,bs);
331:   VecCreate(comm,&xtmp);
332:   VecSetSizes(xtmp,m,PETSC_DECIDE);
333:   VecSetBlockSize(xtmp,bs);
334:   VecSetType(xtmp,vectype);
335:   VecScatterCreate(x,isin,xtmp,NULL,&scatter);

337:   PetscObjectCompose((PetscObject)dmf,"isin",(PetscObject)isin);
338:   PetscObjectCompose((PetscObject)dmf,"scatter",(PetscObject)scatter);
339:   PetscObjectCompose((PetscObject)dmf,"xtmp",(PetscObject)xtmp);
340:   PetscObjectCompose((PetscObject)dmf,"xp",(PetscObject)xp);

342:   VecDestroy(&xred);
343:   VecDestroy(&yred);
344:   VecDestroy(&x);
345:   return 0;
346: }

348: PetscErrorCode DMCreateMatrix_ShellDA(DM dm,Mat *A)
349: {
350:   DM             da;
351:   MPI_Comm       comm;
352:   PetscMPIInt    size;
353:   UserContext    *ctx = NULL;
354:   PetscInt       M,N;

357:   DMShellGetContext(dm,&da);
358:   PetscObjectGetComm((PetscObject)da,&comm);
359:   MPI_Comm_size(comm,&size);
360:   DMCreateMatrix(da,A);
361:   MatGetSize(*A,&M,&N);
362:   PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA (%D x %D)\n",(PetscInt)size,M,N);

364:   DMGetApplicationContext(dm,&ctx);
365:   if (ctx->bcType == NEUMANN) {
366:     MatNullSpace nullspace = NULL;
367:     PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: using neumann bcs\n",(PetscInt)size);

369:     MatGetNullSpace(*A,&nullspace);
370:     if (!nullspace) {
371:       PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: operator does not have nullspace - attaching\n",(PetscInt)size);
372:       MatNullSpaceCreate(comm,PETSC_TRUE,0,0,&nullspace);
373:       MatSetNullSpace(*A,nullspace);
374:       MatNullSpaceDestroy(&nullspace);
375:     } else {
376:       PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: operator already has a nullspace\n",(PetscInt)size);
377:     }
378:   }
379:   return 0;
380: }

382: PetscErrorCode DMCreateGlobalVector_ShellDA(DM dm,Vec *x)
383: {
384:   DM             da;
386:   DMShellGetContext(dm,&da);
387:   DMCreateGlobalVector(da,x);
388:   VecSetDM(*x,dm);
389:   return 0;
390: }

392: PetscErrorCode DMCreateLocalVector_ShellDA(DM dm,Vec *x)
393: {
394:   DM             da;
396:   DMShellGetContext(dm,&da);
397:   DMCreateLocalVector(da,x);
398:   VecSetDM(*x,dm);
399:   return 0;
400: }

402: PetscErrorCode DMCoarsen_ShellDA(DM dm,MPI_Comm comm,DM *dmc)
403: {
405:   *dmc = NULL;
406:   DMGetCoarseDM(dm,dmc);
407:   if (!*dmc) {
408:     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"The coarse DM should never be NULL. The DM hierarchy should have already been defined");
409:   } else {
410:     PetscObjectReference((PetscObject)(*dmc));
411:   }
412:   return 0;
413: }

415: PetscErrorCode DMCreateInterpolation_ShellDA(DM dm1,DM dm2,Mat *mat,Vec *vec)
416: {
417:   DM             da1,da2;
419:   DMShellGetContext(dm1,&da1);
420:   DMShellGetContext(dm2,&da2);
421:   DMCreateInterpolation(da1,da2,mat,vec);
422:   return 0;
423: }

425: PetscErrorCode DMShellDASetUp_TelescopeDMScatter(DM dmf_shell,DM dmc_shell)
426: {
427:   Mat            P = NULL;
428:   DM             dmf = NULL,dmc = NULL;

431:   DMShellGetContext(dmf_shell,&dmf);
432:   if (dmc_shell) {
433:     DMShellGetContext(dmc_shell,&dmc);
434:   }
435:   DMDACreatePermutation_2d(dmc,dmf,&P);
436:   PetscObjectCompose((PetscObject)dmf,"P",(PetscObject)P);
437:   PCTelescopeSetUp_dmda_scatters(dmf,dmc);
438:   PetscObjectComposeFunction((PetscObject)dmf_shell,"PCTelescopeFieldScatter",DMFieldScatter_ShellDA);
439:   PetscObjectComposeFunction((PetscObject)dmf_shell,"PCTelescopeStateScatter",DMStateScatter_ShellDA);
440:   return 0;
441: }

443: PetscErrorCode DMShellDAFieldScatter_Forward(DM dmf,Vec x,DM dmc,Vec xc)
444: {
445:   Mat               P = NULL;
446:   Vec               xp = NULL,xtmp = NULL;
447:   VecScatter        scatter = NULL;
448:   const PetscScalar *x_array;
449:   PetscInt          i,st,ed;

452:   PetscObjectQuery((PetscObject)dmf,"P",(PetscObject*)&P);
453:   PetscObjectQuery((PetscObject)dmf,"xp",(PetscObject*)&xp);
454:   PetscObjectQuery((PetscObject)dmf,"scatter",(PetscObject*)&scatter);
455:   PetscObjectQuery((PetscObject)dmf,"xtmp",(PetscObject*)&xtmp);

461:   MatMultTranspose(P,x,xp);

463:   /* pull in vector x->xtmp */
464:   VecScatterBegin(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);
465:   VecScatterEnd(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);

467:   /* copy vector entries into xred */
468:   VecGetArrayRead(xtmp,&x_array);
469:   if (xc) {
470:     PetscScalar *LA_xred;
471:     VecGetOwnershipRange(xc,&st,&ed);

473:     VecGetArray(xc,&LA_xred);
474:     for (i=0; i<ed-st; i++) {
475:       LA_xred[i] = x_array[i];
476:     }
477:     VecRestoreArray(xc,&LA_xred);
478:   }
479:   VecRestoreArrayRead(xtmp,&x_array);
480:   return 0;
481: }

483: PetscErrorCode DMShellDAFieldScatter_Reverse(DM dmf,Vec y,DM dmc,Vec yc)
484: {
485:   Mat            P = NULL;
486:   Vec            xp = NULL,xtmp = NULL;
487:   VecScatter     scatter = NULL;
488:   PetscScalar    *array;
489:   PetscInt       i,st,ed;

492:   PetscObjectQuery((PetscObject)dmf,"P",(PetscObject*)&P);
493:   PetscObjectQuery((PetscObject)dmf,"xp",(PetscObject*)&xp);
494:   PetscObjectQuery((PetscObject)dmf,"scatter",(PetscObject*)&scatter);
495:   PetscObjectQuery((PetscObject)dmf,"xtmp",(PetscObject*)&xtmp);


502:   /* return vector */
503:   VecGetArray(xtmp,&array);
504:   if (yc) {
505:     const PetscScalar *LA_yred;
506:     VecGetOwnershipRange(yc,&st,&ed);
507:     VecGetArrayRead(yc,&LA_yred);
508:     for (i=0; i<ed-st; i++) {
509:       array[i] = LA_yred[i];
510:     }
511:     VecRestoreArrayRead(yc,&LA_yred);
512:   }
513:   VecRestoreArray(xtmp,&array);
514:   VecScatterBegin(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);
515:   VecScatterEnd(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);
516:   MatMult(P,xp,y);
517:   return 0;
518: }

520: PetscErrorCode DMFieldScatter_ShellDA(DM dmf_shell,Vec x,ScatterMode mode,DM dmc_shell,Vec xc)
521: {
522:   DM             dmf = NULL,dmc = NULL;

525:   DMShellGetContext(dmf_shell,&dmf);
526:   if (dmc_shell) {
527:     DMShellGetContext(dmc_shell,&dmc);
528:   }
529:   if (mode == SCATTER_FORWARD) {
530:     DMShellDAFieldScatter_Forward(dmf,x,dmc,xc);
531:   } else if (mode == SCATTER_REVERSE) {
532:     DMShellDAFieldScatter_Reverse(dmf,x,dmc,xc);
533:   } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell),PETSC_ERR_SUP,"Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
534:   return 0;
535: }

537: PetscErrorCode DMStateScatter_ShellDA(DM dmf_shell,ScatterMode mode,DM dmc_shell)
538: {
539:   PetscMPIInt    size_f = 0,size_c = 0;
541:   MPI_Comm_size(PetscObjectComm((PetscObject)dmf_shell),&size_f);
542:   if (dmc_shell) {
543:     MPI_Comm_size(PetscObjectComm((PetscObject)dmc_shell),&size_c);
544:   }
545:   if (mode == SCATTER_FORWARD) {
546:     PetscPrintf(PetscObjectComm((PetscObject)dmf_shell),"User supplied state scatter (fine [size %d]-> coarse [size %d])\n",(int)size_f,(int)size_c);
547:   } else if (mode == SCATTER_REVERSE) {
548:   } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell),PETSC_ERR_SUP,"Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
549:   return 0;
550: }

552: PetscErrorCode DMShellCreate_ShellDA(DM da,DM *dms)
553: {
555:   if (da) {
556:     DMShellCreate(PetscObjectComm((PetscObject)da),dms);
557:     DMShellSetContext(*dms,da);
558:     DMShellSetCreateGlobalVector(*dms,DMCreateGlobalVector_ShellDA);
559:     DMShellSetCreateLocalVector(*dms,DMCreateLocalVector_ShellDA);
560:     DMShellSetCreateMatrix(*dms,DMCreateMatrix_ShellDA);
561:     DMShellSetCoarsen(*dms,DMCoarsen_ShellDA);
562:     DMShellSetCreateInterpolation(*dms,DMCreateInterpolation_ShellDA);
563:   } else {
564:     *dms = NULL;
565:   }
566:   return 0;
567: }

569: PetscErrorCode DMDestroyShellDMDA(DM *_dm)
570: {
571:   DM             dm,da = NULL;

574:   if (!_dm) return 0;
575:   dm = *_dm;
576:   if (!dm) return 0;

578:   DMShellGetContext(dm,&da);
579:   if (da) {
580:     Vec        vec;
581:     VecScatter scatter = NULL;
582:     IS         is = NULL;
583:     Mat        P = NULL;

585:     PetscObjectQuery((PetscObject)da,"P",(PetscObject*)&P);
586:     MatDestroy(&P);

588:     vec = NULL;
589:     PetscObjectQuery((PetscObject)da,"xp",(PetscObject*)&vec);
590:     VecDestroy(&vec);

592:     PetscObjectQuery((PetscObject)da,"scatter",(PetscObject*)&scatter);
593:     VecScatterDestroy(&scatter);

595:     vec = NULL;
596:     PetscObjectQuery((PetscObject)da,"xtmp",(PetscObject*)&vec);
597:     VecDestroy(&vec);

599:     PetscObjectQuery((PetscObject)da,"isin",(PetscObject*)&is);
600:     ISDestroy(&is);

602:     DMDestroy(&da);
603:   }
604:   DMDestroy(&dm);
605:   *_dm = NULL;
606:   return 0;
607: }

609: PetscErrorCode HierarchyCreate_Basic(DM *dm_f,DM *dm_c,UserContext *ctx)
610: {
611:   DM             dm,dmc,dm_shell,dmc_shell;
612:   PetscMPIInt    rank;

615:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
616:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,17,17,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dm);
617:   DMSetFromOptions(dm);
618:   DMSetUp(dm);
619:   DMDASetUniformCoordinates(dm,0,1,0,1,0,0);
620:   DMDASetFieldName(dm,0,"Pressure");
621:   DMShellCreate_ShellDA(dm,&dm_shell);
622:   DMSetApplicationContext(dm_shell,ctx);

624:   dmc = NULL;
625:   dmc_shell = NULL;
626:   if (rank == 0) {
627:     DMDACreate2d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,17,17,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dmc);
628:     DMSetFromOptions(dmc);
629:     DMSetUp(dmc);
630:     DMDASetUniformCoordinates(dmc,0,1,0,1,0,0);
631:     DMDASetFieldName(dmc,0,"Pressure");
632:     DMShellCreate_ShellDA(dmc,&dmc_shell);
633:     DMSetApplicationContext(dmc_shell,ctx);
634:   }

636:   DMSetCoarseDM(dm_shell,dmc_shell);
637:   DMShellDASetUp_TelescopeDMScatter(dm_shell,dmc_shell);

639:   *dm_f = dm_shell;
640:   *dm_c = dmc_shell;
641:   return 0;
642: }

644: PetscErrorCode HierarchyCreate(PetscInt *_nd,PetscInt *_nref,MPI_Comm **_cl,DM **_dl)
645: {
646:   PetscInt       d,k,ndecomps,ncoarsen,found,nx;
647:   PetscInt       levelrefs,*number;
648:   PetscSubcomm   *pscommlist;
649:   MPI_Comm       *commlist;
650:   DM             *dalist,*dmlist;
651:   PetscBool      set;

654:   ndecomps = 1;
655:   PetscOptionsGetInt(NULL,NULL,"-ndecomps",&ndecomps,NULL);
656:   ncoarsen = ndecomps - 1;

659:   levelrefs = 2;
660:   PetscOptionsGetInt(NULL,NULL,"-level_nrefs",&levelrefs,NULL);
661:   PetscMalloc1(ncoarsen+1,&number);
662:   for (k=0; k<ncoarsen+1; k++) {
663:     number[k] = 2;
664:   }
665:   found = ncoarsen;
666:   set = PETSC_FALSE;
667:   PetscOptionsGetIntArray(NULL,NULL,"-level_comm_red_factor",number,&found,&set);
668:   if (set) {
670:   }

672:   PetscMalloc1(ncoarsen+1,&pscommlist);
673:   for (k=0; k<ncoarsen+1; k++) {
674:     pscommlist[k] = NULL;
675:   }

677:   PetscMalloc1(ndecomps,&commlist);
678:   for (k=0; k<ndecomps; k++) {
679:     commlist[k] = MPI_COMM_NULL;
680:   }
681:   PetscMalloc1(ndecomps*levelrefs,&dalist);
682:   PetscMalloc1(ndecomps*levelrefs,&dmlist);
683:   for (k=0; k<ndecomps*levelrefs; k++) {
684:     dalist[k] = NULL;
685:     dmlist[k] = NULL;
686:   }

688:   CommHierarchyCreate(PETSC_COMM_WORLD,ncoarsen,number,pscommlist);
689:   for (k=0; k<ncoarsen; k++) {
690:     if (pscommlist[k]) {
691:       MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);
692:       if (pscommlist[k]->color == 0) {
693:         PetscCommDuplicate(comm_k,&commlist[k],NULL);
694:       }
695:     }
696:   }
697:   PetscCommDuplicate(PETSC_COMM_WORLD,&commlist[ndecomps-1],NULL);

699:   for (k=0; k<ncoarsen; k++) {
700:     if (pscommlist[k]) {
701:       PetscSubcommDestroy(&pscommlist[k]);
702:     }
703:   }

705:   nx = 17;
706:   PetscOptionsGetInt(NULL,NULL,"-m",&nx,NULL);
707:   for (d=0; d<ndecomps; d++) {
708:     DM   dmroot = NULL;
709:     char name[PETSC_MAX_PATH_LEN];

711:     if (commlist[d] != MPI_COMM_NULL) {
712:       DMDACreate2d(commlist[d],DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,nx,nx,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dmroot);
713:       DMSetUp(dmroot);
714:       DMDASetUniformCoordinates(dmroot,0,1,0,1,0,0);
715:       DMDASetFieldName(dmroot,0,"Pressure");
716:       PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"root-decomp-%D",d);
717:       PetscObjectSetName((PetscObject)dmroot,name);
718:       /*DMView(dmroot,PETSC_VIEWER_STDOUT_(commlist[d]));*/
719:     }

721:     dalist[d*levelrefs + 0] = dmroot;
722:     for (k=1; k<levelrefs; k++) {
723:       DM dmref = NULL;

725:       if (commlist[d] != MPI_COMM_NULL) {
726:         DMRefine(dalist[d*levelrefs + (k-1)],MPI_COMM_NULL,&dmref);
727:         PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"ref%D-decomp-%D",k,d);
728:         PetscObjectSetName((PetscObject)dmref,name);
729:         DMDAGetInfo(dmref,NULL,&nx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
730:         /*DMView(dmref,PETSC_VIEWER_STDOUT_(commlist[d]));*/
731:       }
732:       dalist[d*levelrefs + k] = dmref;
733:     }
734:     MPI_Allreduce(MPI_IN_PLACE,&nx,1,MPIU_INT,MPI_MAX,PETSC_COMM_WORLD);
735:   }

737:   /* create the hierarchy of DMShell's */
738:   for (d=0; d<ndecomps; d++) {
739:     char name[PETSC_MAX_PATH_LEN];

741:     UserContext *ctx = NULL;
742:     if (commlist[d] != MPI_COMM_NULL) {
743:       UserContextCreate(commlist[d],&ctx);
744:       for (k=0; k<levelrefs; k++) {
745:         DMShellCreate_ShellDA(dalist[d*levelrefs + k],&dmlist[d*levelrefs + k]);
746:         DMSetApplicationContext(dmlist[d*levelrefs + k],ctx);
747:         PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"level%D-decomp-%D",k,d);
748:         PetscObjectSetName((PetscObject)dmlist[d*levelrefs + k],name);
749:       }
750:     }
751:   }

753:   /* set all the coarse DMs */
754:   for (k=1; k<ndecomps*levelrefs; k++) { /* skip first DM as it doesn't have a coarse representation */
755:     DM dmfine = NULL,dmcoarse = NULL;

757:     dmfine = dmlist[k];
758:     dmcoarse = dmlist[k-1];
759:     if (dmfine) {
760:       DMSetCoarseDM(dmfine,dmcoarse);
761:     }
762:   }

764:   /* do special setup on the fine DM coupling different decompositions */
765:   for (d=1; d<ndecomps; d++) { /* skip first decomposition as it doesn't have a coarse representation */
766:     DM dmfine = NULL,dmcoarse = NULL;

768:     dmfine = dmlist[d*levelrefs + 0];
769:     dmcoarse = dmlist[(d-1)*levelrefs + (levelrefs-1)];
770:     if (dmfine) {
771:       DMShellDASetUp_TelescopeDMScatter(dmfine,dmcoarse);
772:     }
773:   }

775:   PetscFree(number);
776:   for (k=0; k<ncoarsen; k++) {
777:     PetscSubcommDestroy(&pscommlist[k]);
778:   }
779:   PetscFree(pscommlist);

781:   if (_nd) {
782:     *_nd = ndecomps;
783:   }
784:   if (_nref) {
785:     *_nref = levelrefs;
786:   }
787:   if (_cl) {
788:     *_cl = commlist;
789:   } else {
790:     for (k=0; k<ndecomps; k++) {
791:       if (commlist[k] != MPI_COMM_NULL) {
792:         PetscCommDestroy(&commlist[k]);
793:       }
794:     }
795:     PetscFree(commlist);
796:   }
797:   if (_dl) {
798:     *_dl = dmlist;
799:     PetscFree(dalist);
800:   } else {
801:     for (k=0; k<ndecomps*levelrefs; k++) {
802:       DMDestroy(&dmlist[k]);
803:       DMDestroy(&dalist[k]);
804:     }
805:     PetscFree(dmlist);
806:     PetscFree(dalist);
807:   }
808:   return 0;
809: }

811: PetscErrorCode test_hierarchy(void)
812: {
813:   PetscInt       d,k,nd,nref;
814:   MPI_Comm       *comms;
815:   DM             *dms;

818:   HierarchyCreate(&nd,&nref,&comms,&dms);

820:   /* destroy user context */
821:   for (d=0; d<nd; d++) {
822:     DM first = dms[d*nref+0];

824:     if (first) {
825:       UserContext *ctx = NULL;

827:       DMGetApplicationContext(first,&ctx);
828:       if (ctx) PetscFree(ctx);
829:       DMSetApplicationContext(first,NULL);
830:     }
831:     for (k=1; k<nref; k++) {
832:       DM dm = dms[d*nref+k];
833:       if (dm) {
834:         DMSetApplicationContext(dm,NULL);
835:       }
836:     }
837:   }

839:   /* destroy DMs */
840:   for (k=0; k<nd*nref; k++) {
841:     if (dms[k]) {
842:       DMDestroyShellDMDA(&dms[k]);
843:     }
844:   }
845:   PetscFree(dms);

847:   /* destroy communicators */
848:   for (k=0; k<nd; k++) {
849:     if (comms[k] != MPI_COMM_NULL) {
850:       PetscCommDestroy(&comms[k]);
851:     }
852:   }
853:   PetscFree(comms);
854:   return 0;
855: }

857: PetscErrorCode test_basic(void)
858: {
859:   DM             dmF,dmdaF = NULL,dmC = NULL;
860:   Mat            A;
861:   Vec            x,b;
862:   KSP            ksp;
863:   PC             pc;
864:   UserContext    *user = NULL;

867:   UserContextCreate(PETSC_COMM_WORLD,&user);
868:   HierarchyCreate_Basic(&dmF,&dmC,user);
869:   DMShellGetContext(dmF,&dmdaF);

871:   DMCreateMatrix(dmF,&A);
872:   DMCreateGlobalVector(dmF,&x);
873:   DMCreateGlobalVector(dmF,&b);
874:   ComputeRHS_DMDA(dmdaF,b,user);

876:   KSPCreate(PETSC_COMM_WORLD,&ksp);
877:   KSPSetComputeOperators(ksp,ComputeMatrix_ShellDA,user);
878:   /*KSPSetOperators(ksp,A,A);*/
879:   KSPSetDM(ksp,dmF);
880:   KSPSetDMActive(ksp,PETSC_TRUE);
881:   KSPSetFromOptions(ksp);
882:   KSPGetPC(ksp,&pc);
883:   PCTelescopeSetUseCoarseDM(pc,PETSC_TRUE);

885:   KSPSolve(ksp,b,x);

887:   if (dmC) {
888:     DMDestroyShellDMDA(&dmC);
889:   }
890:   DMDestroyShellDMDA(&dmF);
891:   KSPDestroy(&ksp);
892:   MatDestroy(&A);
893:   VecDestroy(&b);
894:   VecDestroy(&x);
895:   PetscFree(user);
896:   return 0;
897: }

899: PetscErrorCode test_mg(void)
900: {
901:   DM             dmF,dmdaF = NULL,*dms = NULL;
902:   Mat            A;
903:   Vec            x,b;
904:   KSP            ksp;
905:   PetscInt       k,d,nd,nref;
906:   MPI_Comm       *comms = NULL;
907:   UserContext    *user = NULL;

910:   HierarchyCreate(&nd,&nref,&comms,&dms);
911:   dmF = dms[nd*nref-1];

913:   DMShellGetContext(dmF,&dmdaF);
914:   DMGetApplicationContext(dmF,&user);

916:   DMCreateMatrix(dmF,&A);
917:   DMCreateGlobalVector(dmF,&x);
918:   DMCreateGlobalVector(dmF,&b);
919:   ComputeRHS_DMDA(dmdaF,b,user);

921:   KSPCreate(PETSC_COMM_WORLD,&ksp);
922:   KSPSetComputeOperators(ksp,ComputeMatrix_ShellDA,user);
923:   /*KSPSetOperators(ksp,A,A);*/
924:   KSPSetDM(ksp,dmF);
925:   KSPSetDMActive(ksp,PETSC_TRUE);
926:   KSPSetFromOptions(ksp);

928:   KSPSolve(ksp,b,x);

930:   for (d=0; d<nd; d++) {
931:     DM first = dms[d*nref+0];

933:     if (first) {
934:       UserContext *ctx = NULL;

936:       DMGetApplicationContext(first,&ctx);
937:       if (ctx) PetscFree(ctx);
938:       DMSetApplicationContext(first,NULL);
939:     }
940:     for (k=1; k<nref; k++) {
941:       DM dm = dms[d*nref+k];
942:       if (dm) {
943:         DMSetApplicationContext(dm,NULL);
944:       }
945:     }
946:   }

948:   for (k=0; k<nd*nref; k++) {
949:     if (dms[k]) {
950:       DMDestroyShellDMDA(&dms[k]);
951:     }
952:   }
953:   PetscFree(dms);

955:   KSPDestroy(&ksp);
956:   MatDestroy(&A);
957:   VecDestroy(&b);
958:   VecDestroy(&x);

960:   for (k=0; k<nd; k++) {
961:     if (comms[k] != MPI_COMM_NULL) {
962:       PetscCommDestroy(&comms[k]);
963:     }
964:   }
965:   PetscFree(comms);
966:   return 0;
967: }

969: int main(int argc,char **argv)
970: {
971:   PetscInt       test_id = 0;

973:   PetscInitialize(&argc,&argv,(char*)0,help);
974:   PetscOptionsGetInt(NULL,NULL,"-tid",&test_id,NULL);
975:   switch (test_id) {
976:   case 0:
977:     test_basic();
978:       break;
979:   case 1:
980:     test_hierarchy();
981:       break;
982:   case 2:
983:     test_mg();
984:       break;
985:   default:
986:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"-tid must be 0,1,2");
987:   }
988:   PetscFinalize();
989:   return 0;
990: }

992: PetscErrorCode ComputeRHS_DMDA(DM da,Vec b,void *ctx)
993: {
994:   UserContext    *user = (UserContext*)ctx;
995:   PetscInt       i,j,mx,my,xm,ym,xs,ys;
996:   PetscScalar    Hx,Hy;
997:   PetscScalar    **array;
998:   PetscBool      isda = PETSC_FALSE;

1001:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
1003:   DMDAGetInfo(da,NULL,&mx,&my,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1004:   Hx   = 1.0 / (PetscReal)(mx-1);
1005:   Hy   = 1.0 / (PetscReal)(my-1);
1006:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
1007:   DMDAVecGetArray(da,b,&array);
1008:   for (j=ys; j<ys+ym; j++) {
1009:     for (i=xs; i<xs+xm; i++) {
1010:       array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
1011:     }
1012:   }
1013:   DMDAVecRestoreArray(da, b, &array);
1014:   VecAssemblyBegin(b);
1015:   VecAssemblyEnd(b);

1017:   /* force right hand side to be consistent for singular matrix */
1018:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
1019:   if (user->bcType == NEUMANN) {
1020:     MatNullSpace nullspace;

1022:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
1023:     MatNullSpaceRemove(nullspace,b);
1024:     MatNullSpaceDestroy(&nullspace);
1025:   }
1026:   return 0;
1027: }

1029: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
1030: {
1032:   if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
1033:     *rho = centerRho;
1034:   } else {
1035:     *rho = 1.0;
1036:   }
1037:   return 0;
1038: }

1040: PetscErrorCode ComputeMatrix_DMDA(DM da,Mat J,Mat jac,void *ctx)
1041: {
1042:   UserContext    *user = (UserContext*)ctx;
1043:   PetscReal      centerRho;
1044:   PetscInt       i,j,mx,my,xm,ym,xs,ys;
1045:   PetscScalar    v[5];
1046:   PetscReal      Hx,Hy,HydHx,HxdHy,rho;
1047:   MatStencil     row, col[5];
1048:   PetscBool      isda = PETSC_FALSE;

1051:   PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
1053:   MatZeroEntries(jac);
1054:   centerRho = user->rho;
1055:   DMDAGetInfo(da,NULL,&mx,&my,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1056:   Hx        = 1.0 / (PetscReal)(mx-1);
1057:   Hy        = 1.0 / (PetscReal)(my-1);
1058:   HxdHy     = Hx/Hy;
1059:   HydHx     = Hy/Hx;
1060:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
1061:   for (j=ys; j<ys+ym; j++) {
1062:     for (i=xs; i<xs+xm; i++) {
1063:       row.i = i; row.j = j;
1064:       ComputeRho(i, j, mx, my, centerRho, &rho);
1065:       if (i==0 || j==0 || i==mx-1 || j==my-1) {
1066:         if (user->bcType == DIRICHLET) {
1067:           v[0] = 2.0*rho*(HxdHy + HydHx);
1068:           MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
1069:         } else if (user->bcType == NEUMANN) {
1070:           PetscInt numx = 0, numy = 0, num = 0;
1071:           if (j!=0) {
1072:             v[num] = -rho*HxdHy;  col[num].i = i;    col[num].j = j-1;
1073:             numy++; num++;
1074:           }
1075:           if (i!=0) {
1076:             v[num] = -rho*HydHx;  col[num].i = i-1;  col[num].j = j;
1077:             numx++; num++;
1078:           }
1079:           if (i!=mx-1) {
1080:             v[num] = -rho*HydHx;  col[num].i = i+1;  col[num].j = j;
1081:             numx++; num++;
1082:           }
1083:           if (j!=my-1) {
1084:             v[num] = -rho*HxdHy;  col[num].i = i;    col[num].j = j+1;
1085:             numy++; num++;
1086:           }
1087:           v[num] = numx*rho*HydHx + numy*rho*HxdHy;  col[num].i = i; col[num].j = j;
1088:           num++;
1089:           MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
1090:         }
1091:       } else {
1092:         v[0] = -rho*HxdHy;              col[0].i = i;   col[0].j = j-1;
1093:         v[1] = -rho*HydHx;              col[1].i = i-1; col[1].j = j;
1094:         v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
1095:         v[3] = -rho*HydHx;              col[3].i = i+1; col[3].j = j;
1096:         v[4] = -rho*HxdHy;              col[4].i = i;   col[4].j = j+1;
1097:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
1098:       }
1099:     }
1100:   }
1101:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1102:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1103:   return 0;
1104: }

1106: PetscErrorCode ComputeMatrix_ShellDA(KSP ksp,Mat J,Mat jac,void *ctx)
1107: {
1108:   DM             dm,da;
1110:   KSPGetDM(ksp,&dm);
1111:   DMShellGetContext(dm,&da);
1112:   ComputeMatrix_DMDA(da,J,jac,ctx);
1113:   return 0;
1114: }

1116: /*TEST

1118:   test:
1119:     suffix: basic_dirichlet
1120:     nsize: 4
1121:     args: -tid 0 -ksp_monitor_short -pc_type telescope -telescope_ksp_max_it 100000 -telescope_pc_type lu -telescope_ksp_type fgmres -telescope_ksp_monitor_short -ksp_type gcr

1123:   test:
1124:     suffix: basic_neumann
1125:     nsize: 4
1126:     requires: !single
1127:     args: -tid 0 -ksp_monitor_short -pc_type telescope -telescope_ksp_max_it 100000 -telescope_pc_type jacobi -telescope_ksp_type fgmres -telescope_ksp_monitor_short -ksp_type gcr -bc_type neumann

1129:   test:
1130:     suffix: mg_2lv_2mg
1131:     nsize: 6
1132:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 2 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -level_comm_red_factor 6 -mg_coarse_telescope_mg_coarse_pc_type lu

1134:   test:
1135:     suffix: mg_3lv_2mg
1136:     nsize: 4
1137:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 3 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_coarse_pc_type telescope -mg_coarse_telescope_mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_mg_coarse_telescope_pc_type lu -m 5

1139:   test:
1140:     suffix: mg_3lv_2mg_customcommsize
1141:     nsize: 12
1142:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 3 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm  -m 5 -level_comm_red_factor 2,6 -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_coarse_pc_type telescope -mg_coarse_telescope_mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_mg_coarse_telescope_pc_type lu

1144:  TEST*/