Actual source code: galerkin.c
petsc-3.11.4 2019-09-28
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: }