Actual source code: damg.c

  1: #define PETSCSNES_DLL
  2: 
 3:  #include petscda.h
 4:  #include petscksp.h
 5:  #include petscmg.h
 6:  #include petscdmmg.h
 7:  #include private/pcimpl.h

  9: /*
 10:    Code for almost fully managing multigrid/multi-level linear solvers for DA grids
 11: */

 15: /*@C
 16:     DMMGCreate - Creates a DA based multigrid solver object. This allows one to 
 17:       easily implement MG methods on regular grids.

 19:     Collective on MPI_Comm

 21:     Input Parameter:
 22: +   comm - the processors that will share the grids and solution process
 23: .   nlevels - number of multigrid levels 
 24: -   user - an optional user context

 26:     Output Parameters:
 27: .    - the context

 29:     Options Database:
 30: +     -dmmg_nlevels <levels> - number of levels to use
 31: .     -dmmg_galerkin - use Galerkin approach to compute coarser matrices
 32: -     -dmmg_mat_type <type> - matrix type that DMMG should create, defaults to MATAIJ

 34:     Notes:
 35:       To provide a different user context for each level call DMMGSetUser() after calling
 36:       this routine

 38:     Level: advanced

 40: .seealso DMMGDestroy(), DMMGSetUser(), DMMGGetUser(), DMMGSetMatType(), DMMGSetUseGalerkin(), DMMGSetNullSpace(), DMMGSetInitialGuess()

 42: @*/
 43: PetscErrorCode  DMMGCreate(MPI_Comm comm,PetscInt nlevels,void *user,DMMG **dmmg)
 44: {
 46:   PetscInt       i;
 47:   DMMG           *p;
 48:   PetscTruth     galerkin,ftype;
 49:   char           mtype[256];

 52:   PetscOptionsGetInt(0,"-dmmg_nlevels",&nlevels,PETSC_IGNORE);
 53:   PetscOptionsHasName(0,"-dmmg_galerkin",&galerkin);

 55:   PetscMalloc(nlevels*sizeof(DMMG),&p);
 56:   for (i=0; i<nlevels; i++) {
 57:     PetscNew(struct _n_DMMG,&p[i]);
 58:     p[i]->nlevels  = nlevels - i;
 59:     p[i]->comm     = comm;
 60:     p[i]->user     = user;
 61:     p[i]->galerkin = galerkin;
 62:     p[i]->mtype    = MATAIJ;
 63:   }
 64:   p[nlevels-1]->galerkin = PETSC_FALSE;
 65:   *dmmg = p;

 67:   PetscOptionsGetString(PETSC_NULL,"-dmmg_mat_type",mtype,256,&ftype);
 68:   if (ftype) {
 69:     DMMGSetMatType(*dmmg,mtype);
 70:   }
 71:   return(0);
 72: }

 76: /*@C
 77:     DMMGSetMatType - Sets the type of matrices that DMMG will create for its solvers.

 79:     Collective on MPI_Comm 

 81:     Input Parameters:
 82: +    dmmg - the DMMG object created with DMMGCreate()
 83: -    mtype - the matrix type, defaults to MATAIJ

 85:     Level: intermediate

 87: .seealso DMMGDestroy(), DMMGSetUser(), DMMGGetUser(), DMMGCreate(), DMMGSetNullSpace()

 89: @*/
 90: PetscErrorCode  DMMGSetMatType(DMMG *dmmg,MatType mtype)
 91: {
 92:   PetscInt i;
 93: 
 95:   for (i=0; i<dmmg[0]->nlevels; i++) {
 96:     dmmg[i]->mtype  = mtype;
 97:   }
 98:   return(0);
 99: }

103: /*@C
104:     DMMGSetUseGalerkinCoarse - Courses the DMMG to use R*A_f*R^T to form
105:        the coarser matrices from finest 

107:     Collective on DMMG

109:     Input Parameter:
110: .    dmmg - the context

112:     Options Database Keys:
113: .    -dmmg_galerkin 

115:     Level: advanced

117:     Notes: After you have called this you can manually set dmmg[0]->galerkin = PETSC_FALSE
118:        to have the coarsest grid not compute via Galerkin but still have the intermediate
119:        grids computed via Galerkin.

121:        The default behavior of this should be idential to using -pc_mg_galerkin; this offers
122:        more potential flexibility since you can select exactly which levels are done via
123:        Galerkin and which are done via user provided function.

125: .seealso DMMGCreate(), PCMGSetUseGalerkin(), DMMGSetMatType(), DMMGSetNullSpace()

127: @*/
128: PetscErrorCode  DMMGSetUseGalerkinCoarse(DMMG* dmmg)
129: {
130:   PetscInt  i,nlevels = dmmg[0]->nlevels;

133:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");

135:   for (i=0; i<nlevels-1; i++) {
136:     dmmg[i]->galerkin = PETSC_TRUE;
137:   }
138:   return(0);
139: }

143: /*@C
144:     DMMGDestroy - Destroys a DA based multigrid solver object. 

146:     Collective on DMMG

148:     Input Parameter:
149: .    - the context

151:     Level: advanced

153: .seealso DMMGCreate()

155: @*/
156: PetscErrorCode  DMMGDestroy(DMMG *dmmg)
157: {
159:   PetscInt       i,nlevels = dmmg[0]->nlevels;

162:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");

164:   for (i=1; i<nlevels; i++) {
165:     if (dmmg[i]->R) {MatDestroy(dmmg[i]->R);}
166:   }
167:   for (i=0; i<nlevels; i++) {
168:     if (dmmg[i]->dm)      {DMDestroy(dmmg[i]->dm);}
169:     if (dmmg[i]->x)       {VecDestroy(dmmg[i]->x);}
170:     if (dmmg[i]->b)       {VecDestroy(dmmg[i]->b);}
171:     if (dmmg[i]->r)       {VecDestroy(dmmg[i]->r);}
172:     if (dmmg[i]->work1)   {VecDestroy(dmmg[i]->work1);}
173:     if (dmmg[i]->w)       {VecDestroy(dmmg[i]->w);}
174:     if (dmmg[i]->work2)   {VecDestroy(dmmg[i]->work2);}
175:     if (dmmg[i]->lwork1)  {VecDestroy(dmmg[i]->lwork1);}
176:     if (dmmg[i]->B && dmmg[i]->B != dmmg[i]->J) {MatDestroy(dmmg[i]->B);}
177:     if (dmmg[i]->J)         {MatDestroy(dmmg[i]->J);}
178:     if (dmmg[i]->Rscale)    {VecDestroy(dmmg[i]->Rscale);}
179:     if (dmmg[i]->fdcoloring){MatFDColoringDestroy(dmmg[i]->fdcoloring);}
180:     if (dmmg[i]->ksp && !dmmg[i]->snes) {KSPDestroy(dmmg[i]->ksp);}
181:     if (dmmg[i]->snes)      {PetscObjectDestroy((PetscObject)dmmg[i]->snes);}
182:     if (dmmg[i]->inject)    {VecScatterDestroy(dmmg[i]->inject);}
183:     PetscFree(dmmg[i]);
184:   }
185:   PetscFree(dmmg);
186:   return(0);
187: }

191: /*@C
192:     DMMGSetDM - Sets the coarse grid information for the grids

194:     Collective on DMMG

196:     Input Parameter:
197: +   dmmg - the context
198: -   dm - the DA or VecPack object

200:     Options Database Keys:
201: +   -dmmg_refine: Use the input problem as the coarse level and refine. Otherwise, use it as the fine level and coarsen.
202: -   -dmmg_hierarchy: Construct all grids at once

204:     Level: advanced

206: .seealso DMMGCreate(), DMMGDestroy(), DMMGSetUseGalerkin(), DMMGSetMatType()

208: @*/
209: PetscErrorCode  DMMGSetDM(DMMG *dmmg, DM dm)
210: {
211:   PetscInt       nlevels     = dmmg[0]->nlevels;
212:   PetscTruth     doRefine    = PETSC_TRUE;
213:   PetscTruth     doHierarchy = PETSC_FALSE;
214:   PetscInt       i;

218:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");

220:   /* Create DM data structure for all the levels */
221:   PetscOptionsGetTruth(PETSC_NULL, "-dmmg_refine", &doRefine, PETSC_IGNORE);
222:   PetscOptionsHasName(PETSC_NULL, "-dmmg_hierarchy", &doHierarchy);
223:   PetscObjectReference((PetscObject) dm);
224:   if (doRefine) {
225:     dmmg[0]->dm = dm;
226:     if (doHierarchy) {
227: /*       DM *hierarchy; */

229: /*       DMRefineHierarchy(dm, nlevels-1, &hierarchy); */
230: /*       for(i = 1; i < nlevels; ++i) { */
231: /*         dmmg[i]->dm = hierarchy[i-1]; */
232: /*       } */
233:       SETERRQ(PETSC_ERR_SUP, "Refinement hierarchy not yet implemented");
234:     } else {
235:       for(i = 1; i < nlevels; ++i) {
236:         DMRefine(dmmg[i-1]->dm, dmmg[i]->comm, &dmmg[i]->dm);
237:       }
238:     }
239:   } else {
240:     dmmg[nlevels-1]->dm = dm;
241:     if (doHierarchy) {
242:       DM *hierarchy;

244:       DMCoarsenHierarchy(dm, nlevels-1, &hierarchy);
245:       for(i = 0; i < nlevels-1; ++i) {
246:         dmmg[nlevels-2-i]->dm = hierarchy[i];
247:       }
248:     } else {
249: /*       for(i = nlevels-2; i >= 0; --i) { */
250: /*         DMCoarsen(dmmg[i+1]->dm, dmmg[i]->comm, &dmmg[i]->dm); */
251: /*       } */
252:       SETERRQ(PETSC_ERR_SUP, "Sequential coarsening not yet implemented");
253:     }
254:   }
255:   DMMGSetUp(dmmg);
256:   return(0);
257: }

261: /*@C
262:     DMMGSetUp - Prepares the DMMG to solve a system

264:     Collective on DMMG

266:     Input Parameter:
267: .   dmmg - the context

269:     Level: advanced

271: .seealso DMMGCreate(), DMMGDestroy(), DMMG, DMMGSetSNES(), DMMGSetKSP(), DMMGSolve(), DMMGSetMatType()

273: @*/
274: PetscErrorCode  DMMGSetUp(DMMG *dmmg)
275: {
277:   PetscInt       i,nlevels = dmmg[0]->nlevels;


281:   /* Create work vectors and matrix for each level */
282:   for (i=0; i<nlevels; i++) {
283:     DMCreateGlobalVector(dmmg[i]->dm,&dmmg[i]->x);
284:     VecDuplicate(dmmg[i]->x,&dmmg[i]->b);
285:     VecDuplicate(dmmg[i]->x,&dmmg[i]->r);
286:   }

288:   /* Create interpolation/restriction between levels */
289:   for (i=1; i<nlevels; i++) {
290:     DMGetInterpolation(dmmg[i-1]->dm,dmmg[i]->dm,&dmmg[i]->R,PETSC_NULL);
291:   }

293:   return(0);
294: }

298: /*@C
299:     DMMGSolve - Actually solves the (non)linear system defined with the DMMG

301:     Collective on DMMG

303:     Input Parameter:
304: .   dmmg - the context

306:     Level: advanced

308:     Options Database:
309: +   -dmmg_grid_sequence - use grid sequencing to get the initial solution for each level from the previous
310: -   -dmmg_monitor_solution - display the solution at each iteration

312:      Notes: For linear (KSP) problems may be called more than once, uses the same 
313:     matrices but recomputes the right hand side for each new solve. Call DMMGSetKSP()
314:     to generate new matrices.
315:  
316: .seealso DMMGCreate(), DMMGDestroy(), DMMG, DMMGSetSNES(), DMMGSetKSP(), DMMGSetUp(), DMMGSetMatType()

318: @*/
319: PetscErrorCode  DMMGSolve(DMMG *dmmg)
320: {
322:   PetscInt       i,nlevels = dmmg[0]->nlevels;
323:   PetscTruth     gridseq,vecmonitor,flg;

326:   PetscOptionsHasName(0,"-dmmg_grid_sequence",&gridseq);
327:   PetscOptionsHasName(0,"-dmmg_monitor_solution",&vecmonitor);
328:   if (gridseq) {
329:     if (dmmg[0]->initialguess) {
330:       (*dmmg[0]->initialguess)(dmmg[0],dmmg[0]->x);
331:       if (dmmg[0]->ksp && !dmmg[0]->snes) {
332:         KSPSetInitialGuessNonzero(dmmg[0]->ksp,PETSC_TRUE);
333:       }
334:     }
335:     for (i=0; i<nlevels-1; i++) {
336:       (*dmmg[i]->solve)(dmmg,i);
337:       if (vecmonitor) {
338:         VecView(dmmg[i]->x,PETSC_VIEWER_DRAW_(dmmg[i]->comm));
339:       }
340:       MatInterpolate(dmmg[i+1]->R,dmmg[i]->x,dmmg[i+1]->x);
341:       if (dmmg[i+1]->ksp && !dmmg[i+1]->snes) {
342:         KSPSetInitialGuessNonzero(dmmg[i+1]->ksp,PETSC_TRUE);
343:       }
344:     }
345:   } else {
346:     if (dmmg[nlevels-1]->initialguess) {
347:       (*dmmg[nlevels-1]->initialguess)(dmmg[nlevels-1],dmmg[nlevels-1]->x);
348:     }
349:   }

351:   /*VecView(dmmg[nlevels-1]->x,PETSC_VIEWER_DRAW_WORLD);*/

353:   (*DMMGGetFine(dmmg)->solve)(dmmg,nlevels-1);
354:   if (vecmonitor) {
355:      VecView(dmmg[nlevels-1]->x,PETSC_VIEWER_DRAW_(dmmg[nlevels-1]->comm));
356:   }

358:   PetscOptionsHasName(PETSC_NULL,"-dmmg_view",&flg);
359:   if (flg && !PetscPreLoadingOn) {
360:     PetscViewer viewer;
361:     PetscViewerASCIIGetStdout(dmmg[0]->comm,&viewer);
362:     DMMGView(dmmg,viewer);
363:   }
364:   PetscOptionsHasName(PETSC_NULL,"-dmmg_view_binary",&flg);
365:   if (flg && !PetscPreLoadingOn) {
366:     DMMGView(dmmg,PETSC_VIEWER_BINARY_(dmmg[0]->comm));
367:   }
368:   return(0);
369: }

373: PetscErrorCode  DMMGSolveKSP(DMMG *dmmg,PetscInt level)
374: {

378:   if (dmmg[level]->rhs) {
379:     CHKMEMQ;
380:     (*dmmg[level]->rhs)(dmmg[level],dmmg[level]->b);
381:     CHKMEMQ;
382:   }
383:   KSPSolve(dmmg[level]->ksp,dmmg[level]->b,dmmg[level]->x);
384:   return(0);
385: }

387: /*
388:     For each level (of grid sequencing) this sets the interpolation/restriction and 
389:     work vectors needed by the multigrid preconditioner within the KSP 
390:     (for nonlinear problems the KSP inside the SNES) of that level.

392:     Also sets the KSP monitoring on all the levels if requested by user.

394: */
397: PetscErrorCode  DMMGSetUpLevel(DMMG *dmmg,KSP ksp,PetscInt nlevels)
398: {
399:   PetscErrorCode          ierr;
400:   PetscInt                i;
401:   PC                      pc;
402:   PetscTruth              ismg,monitor,ismf,isshell,ismffd;
403:   KSP                     lksp; /* solver internal to the multigrid preconditioner */
404:   MPI_Comm                *comms,comm;
405:   PetscViewerASCIIMonitor ascii;

408:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");

410:   PetscOptionsHasName(PETSC_NULL,"-dmmg_ksp_monitor",&monitor);
411:   if (monitor) {
412:     PetscObjectGetComm((PetscObject)ksp,&comm);
413:     PetscViewerASCIIMonitorCreate(comm,"stdout",1+dmmg[0]->nlevels-nlevels,&ascii);
414:     KSPMonitorSet(ksp,KSPMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
415:   }

417:   /* use fgmres on outer iteration by default */
418:   KSPSetType(ksp,KSPFGMRES);
419:   KSPGetPC(ksp,&pc);
420:   PCSetType(pc,PCMG);
421:   PetscMalloc(nlevels*sizeof(MPI_Comm),&comms);
422:   for (i=0; i<nlevels; i++) {
423:     comms[i] = dmmg[i]->comm;
424:   }
425:   PCMGSetLevels(pc,nlevels,comms);
426:   PetscFree(comms);
427:    PCMGSetType(pc,PC_MG_FULL);

429:   PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
430:   if (ismg) {
431:     /* set solvers for each level */
432:     for (i=0; i<nlevels; i++) {
433:       if (i < nlevels-1) { /* don't set for finest level, they are set in PCApply_MG()*/
434:         PCMGSetX(pc,i,dmmg[i]->x);
435:         PCMGSetRhs(pc,i,dmmg[i]->b);
436:       }
437:       if (i > 0) {
438:         PCMGSetR(pc,i,dmmg[i]->r);
439:       }
440:       if (monitor) {
441:         PCMGGetSmoother(pc,i,&lksp);
442:         PetscObjectGetComm((PetscObject)lksp,&comm);
443:         PetscViewerASCIIMonitorCreate(comm,"stdout",1+dmmg[0]->nlevels-i,&ascii);
444:         KSPMonitorSet(lksp,KSPMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
445:       }
446:       /* If using a matrix free multiply and did not provide an explicit matrix to build
447:          the preconditioner then must use no preconditioner 
448:       */
449:       PetscTypeCompare((PetscObject)dmmg[i]->B,MATSHELL,&isshell);
450:       PetscTypeCompare((PetscObject)dmmg[i]->B,MATDAAD,&ismf);
451:       PetscTypeCompare((PetscObject)dmmg[i]->B,MATMFFD,&ismffd);
452:       if (isshell || ismf || ismffd) {
453:         PC  lpc;
454:         PCMGGetSmoother(pc,i,&lksp);
455:         KSPGetPC(lksp,&lpc);
456:         PCSetType(lpc,PCNONE);
457:       }
458:     }

460:     /* Set interpolation/restriction between levels */
461:     for (i=1; i<nlevels; i++) {
462:       PCMGSetInterpolate(pc,i,dmmg[i]->R);
463:       PCMGSetRestriction(pc,i,dmmg[i]->R);
464:     }
465:   }
466:   return(0);
467: }

471: /*@C
472:     DMMGSetKSP - Sets the linear solver object that will use the grid hierarchy

474:     Collective on DMMG

476:     Input Parameter:
477: +   dmmg - the context
478: .   func - function to compute linear system matrix on each grid level
479: -   rhs - function to compute right hand side on each level (need only work on the finest grid
480:           if you do not use grid sequencing)

482:     Level: advanced

484:     Notes: For linear problems my be called more than once, reevaluates the matrices if it is called more
485:        than once. Call DMMGSolve() directly several times to solve with the same matrix but different 
486:        right hand sides.
487:    
488: .seealso DMMGCreate(), DMMGDestroy, DMMGSetDM(), DMMGSolve(), DMMGSetMatType()

490: @*/
491: PetscErrorCode  DMMGSetKSP(DMMG *dmmg,PetscErrorCode (*rhs)(DMMG,Vec),PetscErrorCode (*func)(DMMG,Mat,Mat))
492: {
494:   PetscInt       i,nlevels = dmmg[0]->nlevels,level;
495:   PetscTruth     galerkin,ismg;
496:   PC             pc;
497:   KSP            lksp;

500:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");
501:   galerkin = dmmg[nlevels - 2 > 0 ? nlevels - 2 : 0]->galerkin;

503:   if (galerkin) {
504:     DMGetMatrix(dmmg[nlevels-1]->dm,dmmg[nlevels-1]->mtype,&dmmg[nlevels-1]->B);
505:     if (!dmmg[nlevels-1]->J) {
506:       dmmg[nlevels-1]->J = dmmg[nlevels-1]->B;
507:     }
508:     (*func)(dmmg[nlevels-1],dmmg[nlevels-1]->J,dmmg[nlevels-1]->B);
509:     for (i=nlevels-2; i>-1; i--) {
510:       if (dmmg[i]->galerkin) {
511:         MatPtAP(dmmg[i+1]->B,dmmg[i+1]->R,MAT_INITIAL_MATRIX,1.0,&dmmg[i]->B);
512:         if (!dmmg[i]->J) {
513:           dmmg[i]->J = dmmg[i]->B;
514:         }
515:       }
516:     }
517:   }

519:   if (!dmmg[0]->ksp) {
520:     /* create solvers for each level if they don't already exist*/
521:     for (i=0; i<nlevels; i++) {

523:       if (!dmmg[i]->B && !dmmg[i]->galerkin) {
524:         DMGetMatrix(dmmg[i]->dm,dmmg[nlevels-1]->mtype,&dmmg[i]->B);
525:       }
526:       if (!dmmg[i]->J) {
527:         dmmg[i]->J = dmmg[i]->B;
528:       }

530:       KSPCreate(dmmg[i]->comm,&dmmg[i]->ksp);
531:       DMMGSetUpLevel(dmmg,dmmg[i]->ksp,i+1);
532:       KSPSetFromOptions(dmmg[i]->ksp);
533:       dmmg[i]->solve = DMMGSolveKSP;
534:       dmmg[i]->rhs   = rhs;
535:     }
536:   }

538:   /* evalute matrix on each level */
539:   for (i=0; i<nlevels; i++) {
540:     if (!dmmg[i]->galerkin) {
541:       (*func)(dmmg[i],dmmg[i]->J,dmmg[i]->B);
542:     }
543:   }

545:   for (i=0; i<nlevels-1; i++) {
546:     KSPSetOptionsPrefix(dmmg[i]->ksp,"dmmg_");
547:   }

549:   for (level=0; level<nlevels; level++) {
550:     KSPSetOperators(dmmg[level]->ksp,dmmg[level]->J,dmmg[level]->B,SAME_NONZERO_PATTERN);
551:     KSPGetPC(dmmg[level]->ksp,&pc);
552:     PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
553:     if (ismg) {
554:       for (i=0; i<=level; i++) {
555:         PCMGGetSmoother(pc,i,&lksp);
556:         KSPSetOperators(lksp,dmmg[i]->J,dmmg[i]->B,SAME_NONZERO_PATTERN);
557:       }
558:     }
559:   }

561:   return(0);
562: }

566: /*@C
567:     DMMGView - prints information on a DA based multi-level preconditioner

569:     Collective on DMMG and PetscViewer

571:     Input Parameter:
572: +   dmmg - the context
573: -   viewer - the viewer

575:     Level: advanced

577: .seealso DMMGCreate(), DMMGDestroy(), DMMGSetMatType(), DMMGSetUseGalerkin()

579: @*/
580: PetscErrorCode  DMMGView(DMMG *dmmg,PetscViewer viewer)
581: {
583:   PetscInt       i,nlevels = dmmg[0]->nlevels;
584:   PetscMPIInt    flag;
585:   MPI_Comm       comm;
586:   PetscTruth     iascii,isbinary;

591:   PetscObjectGetComm((PetscObject)viewer,&comm);
592:   MPI_Comm_compare(comm,dmmg[0]->comm,&flag);
593:   if (flag != MPI_CONGRUENT && flag != MPI_IDENT) {
594:     SETERRQ(PETSC_ERR_ARG_NOTSAMECOMM,"Different communicators in the DMMG and the PetscViewer");
595:   }

597:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
598:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
599:   if (isbinary) {
600:     for (i=0; i<nlevels; i++) {
601:       MatView(dmmg[i]->J,viewer);
602:     }
603:     for (i=1; i<nlevels; i++) {
604:       MatView(dmmg[i]->R,viewer);
605:     }
606:   } else {
607:     if (iascii) {
608:       PetscViewerASCIIPrintf(viewer,"DMMG Object with %D levels\n",nlevels);
609:     }
610:     for (i=0; i<nlevels; i++) {
611:       PetscViewerASCIIPushTab(viewer);
612:       DMView(dmmg[i]->dm,viewer);
613:       PetscViewerASCIIPopTab(viewer);
614:     }
615:     if (iascii) {
616:       PetscViewerASCIIPrintf(viewer,"%s Object on finest level\n",dmmg[nlevels-1]->ksp ? "KSP" : "SNES");
617:       if (dmmg[nlevels-2 > 0 ? nlevels-2 : 0]->galerkin) {
618:         PetscViewerASCIIPrintf(viewer,"Using Galerkin R^T*A*R process to compute coarser matrices\n");
619:       }
620:       PetscViewerASCIIPrintf(viewer,"Using matrix type %s\n",dmmg[nlevels-1]->mtype);
621:     }
622:     if (dmmg[nlevels-1]->ksp) {
623:       KSPView(dmmg[nlevels-1]->ksp,viewer);
624:     } else {
625:       /* use of PetscObjectView() means we do not have to link with libpetscsnes if SNES is not being used */
626:       PetscObjectView((PetscObject)dmmg[nlevels-1]->snes,viewer);
627:     }
628:   }
629:   return(0);
630: }

634: /*@C
635:     DMMGSetNullSpace - Indicates the null space in the linear operator (this is needed by the linear solver)

637:     Collective on DMMG

639:     Input Parameter:
640: +   dmmg - the context
641: .   has_cnst - is the constant vector in the null space
642: .   n - number of null vectors (excluding the possible constant vector)
643: -   func - a function that fills an array of vectors with the null vectors (must be orthonormal), may be PETSC_NULL

645:     Level: advanced

647: .seealso DMMGCreate(), DMMGDestroy, DMMGSetDM(), DMMGSolve(), MatNullSpaceCreate(), KSPSetNullSpace(), DMMGSetUseGalerkin(), DMMGSetMatType()

649: @*/
650: PetscErrorCode  DMMGSetNullSpace(DMMG *dmmg,PetscTruth has_cnst,PetscInt n,PetscErrorCode (*func)(DMMG,Vec[]))
651: {
653:   PetscInt       i,j,nlevels = dmmg[0]->nlevels;
654:   Vec            *nulls = 0;
655:   MatNullSpace   nullsp;
656:   KSP            iksp;
657:   PC             pc,ipc;
658:   PetscTruth     ismg,isred;

661:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");
662:   if (!dmmg[0]->ksp) SETERRQ(PETSC_ERR_ORDER,"Must call AFTER DMMGSetKSP() or DMMGSetSNES()");
663:   if ((n && !func) || (!n && func)) SETERRQ(PETSC_ERR_ARG_INCOMP,"Both n and func() must be set together");
664:   if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Cannot have negative number of vectors in null space n = %D",n)

666:   for (i=0; i<nlevels; i++) {
667:     if (n) {
668:       VecDuplicateVecs(dmmg[i]->b,n,&nulls);
669:       (*func)(dmmg[i],nulls);
670:     }
671:     MatNullSpaceCreate(dmmg[i]->comm,has_cnst,n,nulls,&nullsp);
672:     KSPSetNullSpace(dmmg[i]->ksp,nullsp);
673:     for (j=i; j<nlevels; j++) {
674:       KSPGetPC(dmmg[j]->ksp,&pc);
675:       PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
676:       if (ismg) {
677:         PCMGGetSmoother(pc,i,&iksp);
678:         KSPSetNullSpace(iksp, nullsp);
679:       }
680:     }
681:     MatNullSpaceDestroy(nullsp);
682:     if (n) {
683:       PetscFree(nulls);
684:     }
685:   }
686:   /* make all the coarse grid solvers have LU shift since they are singular */
687:   for (i=0; i<nlevels; i++) {
688:     KSPGetPC(dmmg[i]->ksp,&pc);
689:     PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
690:     if (ismg) {
691:       PCMGGetSmoother(pc,0,&iksp);
692:       KSPGetPC(iksp,&ipc);
693:       PetscTypeCompare((PetscObject)ipc,PCREDUNDANT,&isred);
694:       if (isred) {
695:         PCRedundantGetPC(ipc,&ipc);
696:       }
697:       PCFactorSetShiftPd(ipc,PETSC_TRUE);
698:     }
699:   }
700:   return(0);
701: }

705: /*@C
706:     DMMGInitialGuessCurrent - Use with DMMGSetInitialGuess() to use the current value in the 
707:        solution vector (obtainable with DMMGGetx() as the initial guess)

709:     Collective on DMMG

711:     Input Parameter:
712: +   dmmg - the context
713: -   vec - dummy argument

715:     Level: intermediate

717: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES(), DMMGSetInitialGuess()

719: @*/
720: PetscErrorCode  DMMGInitialGuessCurrent(DMMG dmmg,Vec vec)
721: {
723:   return(0);
724: }

728: /*@C
729:     DMMGSetInitialGuess - Sets the function that computes an initial guess.

731:     Collective on DMMG

733:     Input Parameter:
734: +   dmmg - the context
735: -   guess - the function

737:     Notes: For nonlinear problems, if this is not set, then the current value in the 
738:              solution vector (obtained with DMMGGetX()) is used. Thus is if you doing 'time
739:              stepping' it will use your current solution as the guess for the next timestep.
740:            If grid sequencing is used (via -dmmg_grid_sequence) then the "guess" function
741:              is used only on the coarsest grid.
742:            For linear problems, if this is not set, then 0 is used as an initial guess.
743:              If you would like the linear solver to also (like the nonlinear solver) use
744:              the current solution vector as the initial guess then use DMMGInitialGuessCurrent()
745:              as the function you pass in

747:     Level: intermediate


750: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES(), DMMGInitialGuessCurrent(), DMMGSetGalekin(), DMMGSetMatType(), DMMGSetNullSpace()

752: @*/
753: PetscErrorCode  DMMGSetInitialGuess(DMMG *dmmg,PetscErrorCode (*guess)(DMMG,Vec))
754: {
755:   PetscInt i,nlevels = dmmg[0]->nlevels;

758:   for (i=0; i<nlevels; i++) {
759:     dmmg[i]->initialguess = guess;
760:   }
761:   return(0);
762: }