Actual source code: galerkin.c
petsc-3.10.5 2019-03-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: 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:
167: Either this or PCGalerkinSetInterpolation() or both must be called
169: Level: Intermediate
171: .keywords: PC, set, Galerkin preconditioner
173: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
174: PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
176: @*/
177: PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R)
178: {
183: PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
184: return(0);
185: }
187: /*@
188: PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner
190: Logically Collective on PC
192: Input Parameter:
193: + pc - the preconditioner context
194: - R - the interpolation operator
196: Notes:
197: Either this or PCGalerkinSetRestriction() or both must be called
199: Level: Intermediate
201: .keywords: PC, set, Galerkin preconditioner
203: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
204: PCGalerkinSetRestriction(), PCGalerkinGetKSP()
206: @*/
207: PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P)
208: {
213: PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
214: return(0);
215: }
217: /*@
218: PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
220: Logically Collective
222: Input Parameter:
223: + pc - the preconditioner context
224: . computeAsub - routine that computes the submatrix from the global matrix
225: - ctx - context used by the routine, or NULL
227: Calling sequence of computeAsub:
228: $ computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx);
230: + PC - the Galerkin PC
231: . A - the matrix in the Galerkin PC
232: . Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed
233: . cAp - the submatrix computed by this routine
234: - ctx - optional user-defined function context
236: Level: Intermediate
238: Notes:
239: Instead of providing this routine you can call PCGalerkinGetKSP() and then KSPSetOperators() to provide the submatrix,
240: but that will not work for multiple KSPSolves with different matrices unless you call it for each solve.
242: 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
243: matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.
245: Developer Notes:
246: If the user does not call this routine nor call PCGalerkinGetKSP() and KSPSetOperators() then PCGalerkin could
247: could automatically compute the submatrix via calls to MatGalerkin() or MatRARt()
249: .keywords: PC, get, Galerkin preconditioner, sub preconditioner
251: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
252: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
254: @*/
255: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
256: {
261: PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx));
262: return(0);
263: }
265: /*@
266: PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.
268: Not Collective
270: Input Parameter:
271: . pc - the preconditioner context
273: Output Parameters:
274: . ksp - the KSP object
276: Level: Intermediate
278: Notes:
279: Once you have called this routine you can call KSPSetOperators() on the resulting ksp to provide the operator for the Galerkin problem,
280: an alternative is to use PCGalerkinSetComputeSubmatrix() to provide a routine that computes the submatrix as needed.
282: .keywords: PC, get, Galerkin preconditioner, sub preconditioner
284: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
285: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinSetComputeSubmatrix()
287: @*/
288: PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp)
289: {
295: PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
296: return(0);
297: }
299: static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc)
300: {
301: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
303: const char *prefix;
304: PetscBool flg;
307: KSPGetOptionsPrefix(jac->ksp,&prefix);
308: PetscStrendswith(prefix,"galerkin_",&flg);
309: if (!flg) {
310: PCGetOptionsPrefix(pc,&prefix);
311: KSPSetOptionsPrefix(jac->ksp,prefix);
312: KSPAppendOptionsPrefix(jac->ksp,"galerkin_");
313: }
315: PetscOptionsHead(PetscOptionsObject,"Galerkin options");
316: if (jac->ksp) {
317: KSPSetFromOptions(jac->ksp);
318: }
319: PetscOptionsTail();
320: return(0);
321: }
323: /* -------------------------------------------------------------------------------------------*/
325: /*MC
326: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
328: $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
329: $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)
331: Level: intermediate
333: Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute
334: the operators automatically.
335: Should there be a prefix for the inner KSP.
336: There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP
338: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
339: PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
341: M*/
343: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
344: {
346: PC_Galerkin *jac;
349: PetscNewLog(pc,&jac);
351: pc->ops->apply = PCApply_Galerkin;
352: pc->ops->setup = PCSetUp_Galerkin;
353: pc->ops->reset = PCReset_Galerkin;
354: pc->ops->destroy = PCDestroy_Galerkin;
355: pc->ops->view = PCView_Galerkin;
356: pc->ops->setfromoptions = PCSetFromOptions_Galerkin;
357: pc->ops->applyrichardson = NULL;
359: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
360: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
361: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
363: pc->data = (void*)jac;
365: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
366: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
367: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
368: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin);
369: return(0);
370: }