Actual source code: lsc.c
1: #include <petsc/private/pcimpl.h>
3: typedef struct {
4: PetscBool allocated, commute, scalediag;
5: KSP kspL, kspMass;
6: Vec Avec0, Avec1, Svec0, scale;
7: Mat L;
8: } PC_LSC;
10: static PetscErrorCode PCLSCAllocate_Private(PC pc)
11: {
12: PC_LSC *lsc = (PC_LSC *)pc->data;
13: Mat A;
15: PetscFunctionBegin;
16: if (lsc->allocated) PetscFunctionReturn(PETSC_SUCCESS);
17: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspL));
18: PetscCall(KSPSetNestLevel(lsc->kspL, pc->kspnestlevel));
19: PetscCall(KSPSetErrorIfNotConverged(lsc->kspL, pc->erroriffailure));
20: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspL, (PetscObject)pc, 1));
21: PetscCall(KSPSetType(lsc->kspL, KSPPREONLY));
22: PetscCall(KSPSetOptionsPrefix(lsc->kspL, ((PetscObject)pc)->prefix));
23: PetscCall(KSPAppendOptionsPrefix(lsc->kspL, "lsc_"));
24: PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, NULL, NULL, NULL));
25: PetscCall(MatCreateVecs(A, &lsc->Avec0, &lsc->Avec1));
26: PetscCall(MatCreateVecs(pc->pmat, &lsc->Svec0, NULL));
27: if (lsc->scalediag) PetscCall(VecDuplicate(lsc->Avec0, &lsc->scale));
29: if (lsc->commute) {
30: PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &lsc->kspMass));
31: PetscCall(KSPSetErrorIfNotConverged(lsc->kspMass, pc->erroriffailure));
32: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lsc->kspMass, (PetscObject)pc, 1));
33: PetscCall(KSPSetType(lsc->kspMass, KSPPREONLY));
34: PetscCall(KSPSetOptionsPrefix(lsc->kspMass, ((PetscObject)pc)->prefix));
35: PetscCall(KSPAppendOptionsPrefix(lsc->kspMass, "lsc_mass_"));
36: } else lsc->kspMass = NULL;
38: lsc->allocated = PETSC_TRUE;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode PCSetUp_LSC(PC pc)
43: {
44: PC_LSC *lsc = (PC_LSC *)pc->data;
45: Mat L, Lp, Qscale;
47: PetscFunctionBegin;
48: PetscCall(PCLSCAllocate_Private(pc));
50: /* Query for L operator */
51: PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_L", (PetscObject *)&L));
52: if (!L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_L", (PetscObject *)&L));
53: PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Lp", (PetscObject *)&Lp));
54: if (!Lp) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Lp", (PetscObject *)&Lp));
56: /* Query for mass operator */
57: PetscCall(PetscObjectQuery((PetscObject)pc->pmat, "LSC_Qscale", (PetscObject *)&Qscale));
58: if (!Qscale) PetscCall(PetscObjectQuery((PetscObject)pc->mat, "LSC_Qscale", (PetscObject *)&Qscale));
60: if (lsc->commute) {
61: PetscCheck(L || Lp, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide an L operator for LSC preconditioning when commuting");
62: if (!L && Lp) L = Lp;
63: else if (L && !Lp) Lp = L;
65: PetscCheck(Qscale, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "The user must provide a Qscale matrix for LSC preconditioning when commuting");
66: } else {
67: if (lsc->scale) {
68: if (!Qscale) PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, &Qscale, NULL, NULL, NULL));
69: PetscCall(MatGetDiagonal(Qscale, lsc->scale));
70: PetscCall(VecReciprocal(lsc->scale));
71: }
72: if (!L) {
73: Mat B, C;
74: PetscCall(MatSchurComplementGetSubMatrices(pc->mat, NULL, NULL, &B, &C, NULL));
75: if (lsc->scale) {
76: Mat CAdiaginv;
77: PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &CAdiaginv));
78: PetscCall(MatDiagonalScale(CAdiaginv, NULL, lsc->scale));
79: if (!lsc->L) PetscCall(MatMatMult(CAdiaginv, B, MAT_INITIAL_MATRIX, PETSC_CURRENT, &lsc->L));
80: else PetscCall(MatMatMult(CAdiaginv, B, MAT_REUSE_MATRIX, PETSC_CURRENT, &lsc->L));
81: PetscCall(MatDestroy(&CAdiaginv));
82: } else {
83: if (!lsc->L) {
84: PetscCall(MatProductCreate(C, B, NULL, &lsc->L));
85: PetscCall(MatProductSetType(lsc->L, MATPRODUCT_AB));
86: PetscCall(MatProductSetFromOptions(lsc->L));
87: PetscCall(MatProductSymbolic(lsc->L));
88: }
89: PetscCall(MatProductNumeric(lsc->L));
90: }
91: Lp = L = lsc->L;
92: }
93: }
95: PetscCall(KSPSetOperators(lsc->kspL, L, Lp));
96: PetscCall(KSPSetFromOptions(lsc->kspL));
97: if (lsc->commute) {
98: PetscCall(KSPSetOperators(lsc->kspMass, Qscale, Qscale));
99: PetscCall(KSPSetFromOptions(lsc->kspMass));
100: }
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: static PetscErrorCode PCApply_LSC(PC pc, Vec x, Vec y)
105: {
106: PC_LSC *lsc = (PC_LSC *)pc->data;
107: Mat A, B, C;
109: PetscFunctionBegin;
110: PetscCall(MatSchurComplementGetSubMatrices(pc->mat, &A, NULL, &B, &C, NULL));
111: if (lsc->commute) {
112: PetscCall(KSPSolve(lsc->kspMass, x, lsc->Svec0));
113: PetscCall(KSPCheckSolve(lsc->kspMass, pc, lsc->Svec0));
114: PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
115: PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
116: PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
117: PetscCall(MatMult(A, lsc->Avec1, lsc->Avec0));
118: PetscCall(KSPSolve(lsc->kspL, lsc->Avec0, lsc->Avec1));
119: PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Avec1));
120: PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
121: PetscCall(KSPSolve(lsc->kspMass, lsc->Svec0, y));
122: PetscCall(KSPCheckSolve(lsc->kspMass, pc, y));
123: } else {
124: PetscCall(KSPSolve(lsc->kspL, x, lsc->Svec0));
125: PetscCall(KSPCheckSolve(lsc->kspL, pc, lsc->Svec0));
126: PetscCall(MatMult(B, lsc->Svec0, lsc->Avec0));
127: if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec0, lsc->Avec0, lsc->scale));
128: PetscCall(MatMult(A, lsc->Avec0, lsc->Avec1));
129: if (lsc->scale) PetscCall(VecPointwiseMult(lsc->Avec1, lsc->Avec1, lsc->scale));
130: PetscCall(MatMult(C, lsc->Avec1, lsc->Svec0));
131: PetscCall(KSPSolve(lsc->kspL, lsc->Svec0, y));
132: PetscCall(KSPCheckSolve(lsc->kspL, pc, y));
133: }
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode PCReset_LSC(PC pc)
138: {
139: PC_LSC *lsc = (PC_LSC *)pc->data;
141: PetscFunctionBegin;
142: PetscCall(VecDestroy(&lsc->Avec0));
143: PetscCall(VecDestroy(&lsc->Avec1));
144: PetscCall(VecDestroy(&lsc->Svec0));
145: PetscCall(KSPDestroy(&lsc->kspL));
146: if (lsc->commute) PetscCall(KSPDestroy(&lsc->kspMass));
147: if (lsc->L) PetscCall(MatDestroy(&lsc->L));
148: if (lsc->scale) PetscCall(VecDestroy(&lsc->scale));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode PCDestroy_LSC(PC pc)
153: {
154: PetscFunctionBegin;
155: PetscCall(PCReset_LSC(pc));
156: PetscCall(PetscFree(pc->data));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems *PetscOptionsObject)
161: {
162: PC_LSC *lsc = (PC_LSC *)pc->data;
164: PetscFunctionBegin;
165: PetscOptionsHeadBegin(PetscOptionsObject, "LSC options");
166: {
167: PetscCall(PetscOptionsBool("-pc_lsc_commute", "Whether to commute the LSC preconditioner in the style of Olshanskii", "None", lsc->commute, &lsc->commute, NULL));
168: PetscCall(PetscOptionsBool("-pc_lsc_scale_diag", "Whether to scale BBt products. Will use the inverse of the diagonal of Qscale or A if the former is not provided.", "None", lsc->scalediag, &lsc->scalediag, NULL));
169: PetscCheck(!lsc->scalediag || !lsc->commute, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "Diagonal-based scaling is not used when doing a commuted LSC. Either do not ask for diagonal-based scaling or use non-commuted LSC in the original style of Elman");
170: }
171: PetscOptionsHeadEnd();
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
176: {
177: PC_LSC *jac = (PC_LSC *)pc->data;
178: PetscBool iascii;
180: PetscFunctionBegin;
181: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
182: if (iascii) {
183: PetscCall(PetscViewerASCIIPushTab(viewer));
184: if (jac->kspL) {
185: PetscCall(KSPView(jac->kspL, viewer));
186: } else {
187: PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC KSP object not yet created, hence cannot display"));
188: }
189: if (jac->commute) {
190: if (jac->kspMass) {
191: PetscCall(KSPView(jac->kspMass, viewer));
192: } else {
193: PetscCall(PetscViewerASCIIPrintf(viewer, "PCLSC Mass KSP object not yet created, hence cannot display"));
194: }
195: }
196: PetscCall(PetscViewerASCIIPopTab(viewer));
197: }
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*MC
202: PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators {cite}`elmanhowleshadidshuttleworthtuminaro2006` {cite}`silvester2001efficient`
204: Options Database Key:
205: + -pc_lsc_commute - Whether to commute the LSC preconditioner in the style of Olshanskii
206: - -pc_lsc_scale_diag - Whether to scale $BB^T$ products. Will use the inverse of the diagonal of $Qscale$ or $A$ if the former is not provided
208: Level: intermediate
210: Notes:
211: This preconditioner will normally be used with `PCFIELDSPLIT` to precondition the Schur complement, but
212: it can be used for any Schur complement system. Consider the Schur complement
214: $$
215: S = A11 - A10 A00^{-1} A01
216: $$
218: `PCLSC` currently doesn't do anything with $A11$, so let's assume it is 0. The idea is that a good approximation to
219: $S^{-1}$ is given by
221: $$
222: (A10 A01)^{-1} A10 A00 A01 (A10 A01)^{-1}
223: $$
225: The product $A10 A01$ can be computed for you, but you can provide it (this is
226: usually more efficient anyway). In the case of incompressible flow, $A10 A01$ is a Laplacian; call it $L$. The current
227: interface is to compose $L$ and a matrix from which to construct its preconditioning $Lp$ on the preconditioning matrix.
229: If you had called `KSPSetOperators`(ksp,S,Sp), $S$ should have type `MATSCHURCOMPLEMENT` and $Sp$ can be any type you
230: like (`PCLSC` doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp".
231: For example, you might have setup code like this
233: .vb
234: PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
235: PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
236: .ve
238: And then your Jacobian assembly would look like
240: .vb
241: PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
242: PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
243: if (L) { assembly L }
244: if (Lp) { assemble Lp }
245: .ve
247: With this, you should be able to choose LSC preconditioning, using e.g. the `PCML` algebraic multigrid to solve with L
248: .vb
249: -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
250: .ve
252: Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.
254: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `Block_Preconditioners`, `PCFIELDSPLIT`,
255: `PCFieldSplitGetSubKSP()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`,
256: `MatCreateSchurComplement()`, `MatCreateSchurComplement()`, `MatSchurComplementSetSubMatrices()`, `MatSchurComplementUpdateSubMatrices()`,
257: `MatSchurComplementSetAinvType()`, `MatGetSchurComplement()`
258: M*/
260: PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
261: {
262: PC_LSC *lsc;
264: PetscFunctionBegin;
265: PetscCall(PetscNew(&lsc));
266: pc->data = (void *)lsc;
268: pc->ops->apply = PCApply_LSC;
269: pc->ops->applytranspose = NULL;
270: pc->ops->setup = PCSetUp_LSC;
271: pc->ops->reset = PCReset_LSC;
272: pc->ops->destroy = PCDestroy_LSC;
273: pc->ops->setfromoptions = PCSetFromOptions_LSC;
274: pc->ops->view = PCView_LSC;
275: pc->ops->applyrichardson = NULL;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }