Actual source code: galerkin.c

petsc-3.9.2 2018-05-20
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:   if (jac->P) {
 29:     MatInterpolate(jac->P,jac->x,y);
 30:   } else {
 31:     MatInterpolate(jac->R,jac->x,y);
 32:   }
 33:   return(0);
 34: }

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

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

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

 68:   return(0);
 69: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

160:    Logically Collective on PC

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

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

168:    Level: Intermediate

170: .keywords: PC, set, Galerkin preconditioner

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

175: @*/
176: PetscErrorCode  PCGalerkinSetRestriction(PC pc,Mat R)
177: {

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

186: /*@
187:    PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner

189:    Logically Collective on PC

191:    Input Parameter:
192: +  pc - the preconditioner context
193: -  R - the interpolation operator

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

197:    Level: Intermediate

199: .keywords: PC, set, Galerkin preconditioner

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

204: @*/
205: PetscErrorCode  PCGalerkinSetInterpolation(PC pc,Mat P)
206: {

211:   PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
212:   return(0);
213: }

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

218:    Logically Collective

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

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

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

234:    Level: Intermediate

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

239:           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
240:           matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.

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

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

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

250: @*/
251: PetscErrorCode  PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
252: {

257:   PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx));
258:   return(0);
259: }

261: /*@
262:    PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.

264:    Not Collective

266:    Input Parameter:
267: .  pc - the preconditioner context

269:    Output Parameters:
270: .  ksp - the KSP object

272:    Level: Intermediate

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

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

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

282: @*/
283: PetscErrorCode  PCGalerkinGetKSP(PC pc,KSP *ksp)
284: {

290:   PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
291:   return(0);
292: }

294: static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc)
295: {
296:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
298:   const char     *prefix;
299:   PetscBool      flg;

302:   KSPGetOptionsPrefix(jac->ksp,&prefix);
303:   PetscStrendswith(prefix,"galerkin_",&flg);
304:   if (!flg) {
305:     PCGetOptionsPrefix(pc,&prefix);
306:     KSPSetOptionsPrefix(jac->ksp,prefix);
307:     KSPAppendOptionsPrefix(jac->ksp,"galerkin_");
308:   }

310:   PetscOptionsHead(PetscOptionsObject,"Galerkin options");
311:   if (jac->ksp) {
312:     KSPSetFromOptions(jac->ksp);
313:   }
314:   PetscOptionsTail();
315:   return(0);
316: }

318: /* -------------------------------------------------------------------------------------------*/

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

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

326:    Level: intermediate

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

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

336: M*/

338: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
339: {
341:   PC_Galerkin    *jac;

344:   PetscNewLog(pc,&jac);

346:   pc->ops->apply           = PCApply_Galerkin;
347:   pc->ops->setup           = PCSetUp_Galerkin;
348:   pc->ops->reset           = PCReset_Galerkin;
349:   pc->ops->destroy         = PCDestroy_Galerkin;
350:   pc->ops->view            = PCView_Galerkin;
351:   pc->ops->setfromoptions  = PCSetFromOptions_Galerkin;
352:   pc->ops->applyrichardson = NULL;

354:   KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
355:   KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
356:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);

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

360:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
361:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
362:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
363:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin);
364:   return(0);
365: }