Actual source code: galerkin.c


  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: {
 18:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

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

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

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

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

 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

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

 74:   PetscFunctionBegin;
 75:   PetscCall(MatDestroy(&jac->R));
 76:   PetscCall(MatDestroy(&jac->P));
 77:   PetscCall(VecDestroy(&jac->x));
 78:   PetscCall(VecDestroy(&jac->b));
 79:   PetscCall(KSPReset(jac->ksp));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

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

 87:   PetscFunctionBegin;
 88:   PetscCall(PCReset_Galerkin(pc));
 89:   PetscCall(KSPDestroy(&jac->ksp));
 90:   PetscCall(PetscFree(pc->data));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
 95: {
 96:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;
 97:   PetscBool    iascii;

 99:   PetscFunctionBegin;
100:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
101:   if (iascii) {
102:     PetscCall(PetscViewerASCIIPrintf(viewer, "  KSP on Galerkin follow\n"));
103:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
104:   }
105:   PetscCall(KSPView(jac->ksp, viewer));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
110: {
111:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

113:   PetscFunctionBegin;
114:   *ksp = jac->ksp;
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
119: {
120:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

122:   PetscFunctionBegin;
123:   PetscCall(PetscObjectReference((PetscObject)R));
124:   PetscCall(MatDestroy(&jac->R));
125:   jac->R = R;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
130: {
131:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

133:   PetscFunctionBegin;
134:   PetscCall(PetscObjectReference((PetscObject)P));
135:   PetscCall(MatDestroy(&jac->P));
136:   jac->P = P;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
141: {
142:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;

144:   PetscFunctionBegin;
145:   jac->computeasub     = computeAsub;
146:   jac->computeasub_ctx = ctx;
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: /*@
151:    PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner

153:    Logically Collective

155:    Input Parameters:
156: +  pc - the preconditioner context
157: -  R - the restriction operator

159:    Level: Intermediate

161:    Note:
162:    Either this or `PCGalerkinSetInterpolation()` or both must be called

164: .seealso: `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
165:           `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
166: @*/
167: PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
168: {
169:   PetscFunctionBegin;
171:   PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: /*@
176:    PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner

178:    Logically Collective

180:    Input Parameters:
181: +  pc - the preconditioner context
182: -  R - the interpolation operator

184:    Level: Intermediate

186:    Note:
187:    Either this or `PCGalerkinSetRestriction()` or both must be called

189: .seealso: `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
190:           `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
191: @*/
192: PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
193: {
194:   PetscFunctionBegin;
196:   PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

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

203:    Logically Collective

205:    Input Parameters:
206: +  pc - the preconditioner context
207: .  computeAsub - routine that computes the submatrix from the global matrix
208: -  ctx - context used by the routine, or `NULL`

210:    Calling sequence of `computeAsub`:
211: $  PetscErrorCode computeAsub(PC pc, Mat A, Mat Ap, Mat *cAP, void *ctx);
212: +  PC - the `PCGALERKIN`
213: .  A - the matrix in the `PCGALERKIN`
214: .  Ap - the computed submatrix from any previous computation, if `NULL` it has not previously been computed
215: .  cAp - the submatrix computed by this routine
216: -  ctx - optional user-defined function context

218:    Level: Intermediate

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

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

227:    Developer Note:
228:    If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN`
229:    could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()`

231: .seealso: `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
232:           `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
233: @*/
234: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
235: {
236:   PetscFunctionBegin;
238:   PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode(*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: /*@
243:    PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`

245:    Not Collective

247:    Input Parameter:
248: .  pc - the preconditioner context

250:    Output Parameter:
251: .  ksp - the `KSP` object

253:    Level: Intermediate

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

259: .seealso: `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
260:           `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
261: @*/
262: PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
263: {
264:   PetscFunctionBegin;
267:   PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
268:   PetscFunctionReturn(PETSC_SUCCESS);
269: }

271: static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems *PetscOptionsObject)
272: {
273:   PC_Galerkin *jac = (PC_Galerkin *)pc->data;
274:   const char  *prefix;
275:   PetscBool    flg;

277:   PetscFunctionBegin;
278:   PetscCall(KSPGetOptionsPrefix(jac->ksp, &prefix));
279:   PetscCall(PetscStrendswith(prefix, "galerkin_", &flg));
280:   if (!flg) {
281:     PetscCall(PCGetOptionsPrefix(pc, &prefix));
282:     PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
283:     PetscCall(KSPAppendOptionsPrefix(jac->ksp, "galerkin_"));
284:   }

286:   PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
287:   if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
288:   PetscOptionsHeadEnd();
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

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

295:     Level: intermediate

297:     Note:
298:     Use
299: .vb
300:      `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P)
301:      `PCGalerkinGetKSP`(pc,&ksp);
302:      `KSPSetOperators`(ksp,A,....)
303:      ...
304: .ve

306:     Developer Notes:
307:     If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute
308:     the operators automatically.

310:     Should there be a prefix for the inner `KSP`?

312:     There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`

314: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
315:           `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
316: M*/

318: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
319: {
320:   PC_Galerkin *jac;

322:   PetscFunctionBegin;
323:   PetscCall(PetscNew(&jac));

325:   pc->ops->apply           = PCApply_Galerkin;
326:   pc->ops->setup           = PCSetUp_Galerkin;
327:   pc->ops->reset           = PCReset_Galerkin;
328:   pc->ops->destroy         = PCDestroy_Galerkin;
329:   pc->ops->view            = PCView_Galerkin;
330:   pc->ops->setfromoptions  = PCSetFromOptions_Galerkin;
331:   pc->ops->applyrichardson = NULL;

333:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
334:   PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
335:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));

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

339:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin));
340:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin));
341:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin));
342:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin));
343:   PetscFunctionReturn(PETSC_SUCCESS);
344: }