Actual source code: galerkin.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  2: /*
  3:       Defines a preconditioner defined by R^T S R
  4: */
  5:  #include <petsc/private/pcimpl.h>
  6:  #include <petscksp.h>

  8: typedef struct {
  9:   KSP            ksp;
 10:   Mat            R,P;
 11:   Vec            b,x;
 12:   PetscErrorCode (*computeasub)(PC,Mat,Mat,Mat*,void*);
 13:   void           *computeasub_ctx;
 14: } PC_Galerkin;

 16: static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y)
 17: {
 19:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

 22:   if (jac->R) {
 23:     MatRestrict(jac->R,x,jac->b);
 24:   } else {
 25:     MatRestrict(jac->P,x,jac->b);
 26:   }
 27:   KSPSolve(jac->ksp,jac->b,jac->x);
 28:   KSPCheckSolve(jac->ksp,pc,jac->x);
 29:   if (jac->P) {
 30:     MatInterpolate(jac->P,jac->x,y);
 31:   } else {
 32:     MatInterpolate(jac->R,jac->x,y);
 33:   }
 34:   return(0);
 35: }

 37: static PetscErrorCode PCSetUp_Galerkin(PC pc)
 38: {
 40:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
 41:   PetscBool      a;
 42:   Vec            *xx,*yy;

 45:   if (jac->computeasub) {
 46:     Mat Ap;
 47:     if (!pc->setupcalled) {
 48:       (*jac->computeasub)(pc,pc->pmat,NULL,&Ap,jac->computeasub_ctx);
 49:       KSPSetOperators(jac->ksp,Ap,Ap);
 50:       MatDestroy(&Ap);
 51:     } else {
 52:       KSPGetOperators(jac->ksp,NULL,&Ap);
 53:       (*jac->computeasub)(pc,pc->pmat,Ap,NULL,jac->computeasub_ctx);
 54:     }
 55:   }

 57:   if (!jac->x) {
 58:     KSPGetOperatorsSet(jac->ksp,&a,NULL);
 59:     if (!a) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
 60:     KSPCreateVecs(jac->ksp,1,&xx,1,&yy);
 61:     jac->x = *xx;
 62:     jac->b = *yy;
 63:     PetscFree(xx);
 64:     PetscFree(yy);
 65:   }
 66:   if (!jac->R && !jac->P) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()");
 67:   /* should check here that sizes of R/P match size of a */

 69:   return(0);
 70: }

 72: static PetscErrorCode PCReset_Galerkin(PC pc)
 73: {
 74:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

 78:   MatDestroy(&jac->R);
 79:   MatDestroy(&jac->P);
 80:   VecDestroy(&jac->x);
 81:   VecDestroy(&jac->b);
 82:   KSPReset(jac->ksp);
 83:   return(0);
 84: }

 86: static PetscErrorCode PCDestroy_Galerkin(PC pc)
 87: {
 88:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

 92:   PCReset_Galerkin(pc);
 93:   KSPDestroy(&jac->ksp);
 94:   PetscFree(pc->data);
 95:   return(0);
 96: }

 98: static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer)
 99: {
100:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
102:   PetscBool      iascii;

105:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
106:   if (iascii) {
107:     PetscViewerASCIIPrintf(viewer,"  KSP on Galerkin follow\n");
108:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
109:   }
110:   KSPView(jac->ksp,viewer);
111:   return(0);
112: }

114: static PetscErrorCode  PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp)
115: {
116:   PC_Galerkin *jac = (PC_Galerkin*)pc->data;

119:   *ksp = jac->ksp;
120:   return(0);
121: }

123: static PetscErrorCode  PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)
124: {
125:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

129:   PetscObjectReference((PetscObject)R);
130:   MatDestroy(&jac->R);
131:   jac->R = R;
132:   return(0);
133: }

135: static PetscErrorCode  PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)
136: {
137:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

141:   PetscObjectReference((PetscObject)P);
142:   MatDestroy(&jac->P);
143:   jac->P = P;
144:   return(0);
145: }

147: static PetscErrorCode  PCGalerkinSetComputeSubmatrix_Galerkin(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
148: {
149:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

152:   jac->computeasub     = computeAsub;
153:   jac->computeasub_ctx = ctx;
154:   return(0);
155: }

157: /* -------------------------------------------------------------------------------- */
158: /*@
159:    PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner

161:    Logically Collective on PC

163:    Input Parameter:
164: +  pc - the preconditioner context
165: -  R - the restriction operator

167:    Notes:
168:     Either this or PCGalerkinSetInterpolation() or both must be called

170:    Level: Intermediate

172: .keywords: PC, set, Galerkin preconditioner

174: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
175:            PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

177: @*/
178: PetscErrorCode  PCGalerkinSetRestriction(PC pc,Mat R)
179: {

184:   PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
185:   return(0);
186: }

188: /*@
189:    PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner

191:    Logically Collective on PC

193:    Input Parameter:
194: +  pc - the preconditioner context
195: -  R - the interpolation operator

197:    Notes:
198:     Either this or PCGalerkinSetRestriction() or both must be called

200:    Level: Intermediate

202: .keywords: PC, set, Galerkin preconditioner

204: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
205:            PCGalerkinSetRestriction(), PCGalerkinGetKSP()

207: @*/
208: PetscErrorCode  PCGalerkinSetInterpolation(PC pc,Mat P)
209: {

214:   PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
215:   return(0);
216: }

218: /*@
219:    PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix

221:    Logically Collective

223:    Input Parameter:
224: +  pc - the preconditioner context
225: .  computeAsub - routine that computes the submatrix from the global matrix
226: -  ctx - context used by the routine, or NULL

228:    Calling sequence of computeAsub:
229: $    computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx);

231: +  PC - the Galerkin PC
232: .  A - the matrix in the Galerkin PC
233: .  Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed
234: .  cAp - the submatrix computed by this routine
235: -  ctx - optional user-defined function context

237:    Level: Intermediate

239:    Notes:
240:     Instead of providing this routine you can call PCGalerkinGetKSP() and then KSPSetOperators() to provide the submatrix,
241:           but that will not work for multiple KSPSolves with different matrices unless you call it for each solve.

243:           This routine is called each time the outter matrix is changed. In the first call the Ap argument is NULL and the routine should create the
244:           matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.

246:    Developer Notes:
247:     If the user does not call this routine nor call PCGalerkinGetKSP() and KSPSetOperators() then PCGalerkin could
248:                     could automatically compute the submatrix via calls to MatGalerkin() or MatRARt()

250: .keywords: PC, get, Galerkin preconditioner, sub preconditioner

252: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
253:            PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

255: @*/
256: PetscErrorCode  PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
257: {

262:   PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx));
263:   return(0);
264: }

266: /*@
267:    PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.

269:    Not Collective

271:    Input Parameter:
272: .  pc - the preconditioner context

274:    Output Parameters:
275: .  ksp - the KSP object

277:    Level: Intermediate

279:    Notes:
280:     Once you have called this routine you can call KSPSetOperators() on the resulting ksp to provide the operator for the Galerkin problem,
281:           an alternative is to use PCGalerkinSetComputeSubmatrix() to provide a routine that computes the submatrix as needed.

283: .keywords: PC, get, Galerkin preconditioner, sub preconditioner

285: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
286:            PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinSetComputeSubmatrix()

288: @*/
289: PetscErrorCode  PCGalerkinGetKSP(PC pc,KSP *ksp)
290: {

296:   PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
297:   return(0);
298: }

300: static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc)
301: {
302:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
304:   const char     *prefix;
305:   PetscBool      flg;

308:   KSPGetOptionsPrefix(jac->ksp,&prefix);
309:   PetscStrendswith(prefix,"galerkin_",&flg);
310:   if (!flg) {
311:     PCGetOptionsPrefix(pc,&prefix);
312:     KSPSetOptionsPrefix(jac->ksp,prefix);
313:     KSPAppendOptionsPrefix(jac->ksp,"galerkin_");
314:   }

316:   PetscOptionsHead(PetscOptionsObject,"Galerkin options");
317:   if (jac->ksp) {
318:     KSPSetFromOptions(jac->ksp);
319:   }
320:   PetscOptionsTail();
321:   return(0);
322: }

324: /* -------------------------------------------------------------------------------------------*/

326: /*MC
327:      PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)

329: $   Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
330: $   PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)

332:    Level: intermediate

334:    Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute
335:                    the operators automatically.
336:                    Should there be a prefix for the inner KSP.
337:                    There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP

339: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
340:            PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

342: M*/

344: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
345: {
347:   PC_Galerkin    *jac;

350:   PetscNewLog(pc,&jac);

352:   pc->ops->apply           = PCApply_Galerkin;
353:   pc->ops->setup           = PCSetUp_Galerkin;
354:   pc->ops->reset           = PCReset_Galerkin;
355:   pc->ops->destroy         = PCDestroy_Galerkin;
356:   pc->ops->view            = PCView_Galerkin;
357:   pc->ops->setfromoptions  = PCSetFromOptions_Galerkin;
358:   pc->ops->applyrichardson = NULL;

360:   KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
361:   KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
362:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);

364:   pc->data = (void*)jac;

366:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
367:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
368:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
369:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin);
370:   return(0);
371: }