Actual source code: galerkin.c
1: /*
2: Defines a preconditioner defined by R^T S R
3: */
4: #include <petsc/private/pcimpl.h>
5: #include <petscksp.h>
7: typedef struct {
8: KSP ksp;
9: Mat R, P;
10: Vec b, x;
11: PetscErrorCode (*computeasub)(PC, Mat, Mat, Mat *, void *);
12: void *computeasub_ctx;
13: } PC_Galerkin;
15: static PetscErrorCode PCApply_Galerkin(PC pc, Vec x, Vec y)
16: {
17: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
19: PetscFunctionBegin;
20: if (jac->R) {
21: PetscCall(MatRestrict(jac->R, x, jac->b));
22: } else {
23: PetscCall(MatRestrict(jac->P, x, jac->b));
24: }
25: PetscCall(KSPSolve(jac->ksp, jac->b, jac->x));
26: PetscCall(KSPCheckSolve(jac->ksp, pc, jac->x));
27: if (jac->P) {
28: PetscCall(MatInterpolate(jac->P, jac->x, y));
29: } else {
30: PetscCall(MatInterpolate(jac->R, jac->x, y));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode PCSetUp_Galerkin(PC pc)
36: {
37: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
38: PetscBool a;
39: Vec *xx, *yy;
41: PetscFunctionBegin;
42: if (jac->computeasub) {
43: Mat Ap;
44: if (!pc->setupcalled) {
45: PetscCall((*jac->computeasub)(pc, pc->pmat, NULL, &Ap, jac->computeasub_ctx));
46: PetscCall(KSPSetOperators(jac->ksp, Ap, Ap));
47: PetscCall(MatDestroy(&Ap));
48: } else {
49: PetscCall(KSPGetOperators(jac->ksp, NULL, &Ap));
50: PetscCall((*jac->computeasub)(pc, pc->pmat, Ap, NULL, jac->computeasub_ctx));
51: }
52: }
54: if (!jac->x) {
55: PetscCall(KSPGetOperatorsSet(jac->ksp, &a, NULL));
56: PetscCheck(a, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
57: PetscCall(KSPCreateVecs(jac->ksp, 1, &xx, 1, &yy));
58: jac->x = *xx;
59: jac->b = *yy;
60: PetscCall(PetscFree(xx));
61: PetscCall(PetscFree(yy));
62: }
63: PetscCheck(jac->R || jac->P, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()");
64: /* should check here that sizes of R/P match size of a */
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode PCReset_Galerkin(PC pc)
70: {
71: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
73: PetscFunctionBegin;
74: PetscCall(MatDestroy(&jac->R));
75: PetscCall(MatDestroy(&jac->P));
76: PetscCall(VecDestroy(&jac->x));
77: PetscCall(VecDestroy(&jac->b));
78: PetscCall(KSPReset(jac->ksp));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode PCDestroy_Galerkin(PC pc)
83: {
84: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
86: PetscFunctionBegin;
87: PetscCall(PCReset_Galerkin(pc));
88: PetscCall(KSPDestroy(&jac->ksp));
89: PetscCall(PetscFree(pc->data));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
94: {
95: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
96: PetscBool iascii;
98: PetscFunctionBegin;
99: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
100: if (iascii) {
101: PetscCall(PetscViewerASCIIPrintf(viewer, " KSP on Galerkin follow\n"));
102: PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
103: }
104: PetscCall(KSPView(jac->ksp, viewer));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
109: {
110: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
112: PetscFunctionBegin;
113: *ksp = jac->ksp;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
118: {
119: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
121: PetscFunctionBegin;
122: PetscCall(PetscObjectReference((PetscObject)R));
123: PetscCall(MatDestroy(&jac->R));
124: jac->R = R;
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
129: {
130: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
132: PetscFunctionBegin;
133: PetscCall(PetscObjectReference((PetscObject)P));
134: PetscCall(MatDestroy(&jac->P));
135: jac->P = P;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
140: {
141: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
143: PetscFunctionBegin;
144: jac->computeasub = computeAsub;
145: jac->computeasub_ctx = ctx;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner
152: Logically Collective
154: Input Parameters:
155: + pc - the preconditioner context
156: - R - the restriction operator
158: Level: intermediate
160: Note:
161: Either this or `PCGalerkinSetInterpolation()` or both must be called
163: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
164: `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
165: @*/
166: PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
167: {
168: PetscFunctionBegin;
170: PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner
177: Logically Collective
179: Input Parameters:
180: + pc - the preconditioner context
181: - P - the interpolation operator
183: Level: intermediate
185: Note:
186: Either this or `PCGalerkinSetRestriction()` or both must be called
188: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
189: `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
190: @*/
191: PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
192: {
193: PetscFunctionBegin;
195: PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@C
200: PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
202: Logically Collective
204: Input Parameters:
205: + pc - the preconditioner context
206: . computeAsub - routine that computes the submatrix from the global matrix
207: - ctx - context used by the routine, or `NULL`
209: Calling sequence of `computeAsub`:
210: + pc - the `PCGALERKIN` preconditioner
211: . A - the matrix in the `PCGALERKIN`
212: . Ap - the computed submatrix from any previous computation, if `NULL` it has not previously been computed
213: . cAp - the submatrix computed by this routine
214: - ctx - optional user-defined function context
216: Level: intermediate
218: Notes:
219: Instead of providing this routine you can call `PCGalerkinGetKSP()` and then `KSPSetOperators()` to provide the submatrix,
220: but that will not work for multiple `KSPSolve()`s with different matrices unless you call it for each solve.
222: 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
223: matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.
225: Developer Notes:
226: If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN`
227: could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()`
229: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
230: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
231: @*/
232: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC pc, Mat A, Mat Ap, Mat *cAp, void *ctx), void *ctx)
233: {
234: PetscFunctionBegin;
236: PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode(*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`
243: Not Collective
245: Input Parameter:
246: . pc - the preconditioner context
248: Output Parameter:
249: . ksp - the `KSP` object
251: Level: intermediate
253: Note:
254: Once you have called this routine you can call `KSPSetOperators()` on the resulting `KSP` to provide the operator for the Galerkin problem,
255: an alternative is to use `PCGalerkinSetComputeSubmatrix()` to provide a routine that computes the submatrix as needed.
257: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
258: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
259: @*/
260: PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
261: {
262: PetscFunctionBegin;
264: PetscAssertPointer(ksp, 2);
265: PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems *PetscOptionsObject)
270: {
271: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
272: const char *prefix;
273: PetscBool flg;
275: PetscFunctionBegin;
276: PetscCall(KSPGetOptionsPrefix(jac->ksp, &prefix));
277: PetscCall(PetscStrendswith(prefix, "galerkin_", &flg));
278: if (!flg) {
279: PetscCall(PCGetOptionsPrefix(pc, &prefix));
280: PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
281: PetscCall(KSPAppendOptionsPrefix(jac->ksp, "galerkin_"));
282: }
284: PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
285: if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
286: PetscOptionsHeadEnd();
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*MC
291: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
293: Level: intermediate
295: Note:
296: Use
297: .vb
298: `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P)
299: `PCGalerkinGetKSP`(pc,&ksp);
300: `KSPSetOperators`(ksp,A,....)
301: ...
302: .ve
304: Developer Notes:
305: If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute
306: the operators automatically.
308: Should there be a prefix for the inner `KSP`?
310: There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`
312: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
313: `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
314: M*/
316: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
317: {
318: PC_Galerkin *jac;
320: PetscFunctionBegin;
321: PetscCall(PetscNew(&jac));
323: pc->ops->apply = PCApply_Galerkin;
324: pc->ops->setup = PCSetUp_Galerkin;
325: pc->ops->reset = PCReset_Galerkin;
326: pc->ops->destroy = PCDestroy_Galerkin;
327: pc->ops->view = PCView_Galerkin;
328: pc->ops->setfromoptions = PCSetFromOptions_Galerkin;
329: pc->ops->applyrichardson = NULL;
331: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
332: PetscCall(KSPSetNestLevel(jac->ksp, pc->kspnestlevel));
333: PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
334: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));
336: pc->data = (void *)jac;
338: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin));
339: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin));
340: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin));
341: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }