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: }