1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
2
3 typedef struct {
4 PetscBool allocated, commute, scalediag;
5 KSP kspL, kspMass;
6 Vec Avec0, Avec1, Svec0, scale;
7 Mat L;
8 } PC_LSC;
9
PCLSCAllocate_Private(PC pc)10 static PetscErrorCode PCLSCAllocate_Private(PC pc)
11 {
12 PC_LSC *lsc = (PC_LSC *)pc->data;
13 Mat A;
14
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));
28
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;
37
38 lsc->allocated = PETSC_TRUE;
39 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41
PCSetUp_LSC(PC pc)42 static PetscErrorCode PCSetUp_LSC(PC pc)
43 {
44 PC_LSC *lsc = (PC_LSC *)pc->data;
45 Mat L, Lp, Qscale;
46
47 PetscFunctionBegin;
48 PetscCall(PCLSCAllocate_Private(pc));
49
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));
55
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));
59
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;
64
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 }
94
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 }
103
PCApply_LSC(PC pc,Vec x,Vec y)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;
108
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 }
136
PCReset_LSC(PC pc)137 static PetscErrorCode PCReset_LSC(PC pc)
138 {
139 PC_LSC *lsc = (PC_LSC *)pc->data;
140
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 }
151
PCDestroy_LSC(PC pc)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 }
159
PCSetFromOptions_LSC(PC pc,PetscOptionItems PetscOptionsObject)160 static PetscErrorCode PCSetFromOptions_LSC(PC pc, PetscOptionItems PetscOptionsObject)
161 {
162 PC_LSC *lsc = (PC_LSC *)pc->data;
163
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 }
174
PCView_LSC(PC pc,PetscViewer viewer)175 static PetscErrorCode PCView_LSC(PC pc, PetscViewer viewer)
176 {
177 PC_LSC *jac = (PC_LSC *)pc->data;
178 PetscBool isascii;
179
180 PetscFunctionBegin;
181 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
182 if (isascii) {
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 }
200
201 /*MC
202 PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators {cite}`elmanhowleshadidshuttleworthtuminaro2006` {cite}`silvester2001efficient`
203
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
207
208 Level: intermediate
209
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
213
214 $$
215 S = A11 - A10 A00^{-1} A01
216 $$
217
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
220
221 $$
222 (A10 A01)^{-1} A10 A00 A01 (A10 A01)^{-1}
223 $$
224
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 a preconditioner $Lp$ on the matrix.
228
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
232
233 .vb
234 PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L);
235 PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp);
236 .ve
237
238 And then your Jacobian assembly would look like
239
240 .vb
241 PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L);
242 PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp);
243 if (L) { assemble L }
244 if (Lp) { assemble Lp }
245 .ve
246
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
251
252 Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners.
253
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*/
259
PCCreate_LSC(PC pc)260 PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc)
261 {
262 PC_LSC *lsc;
263
264 PetscFunctionBegin;
265 PetscCall(PetscNew(&lsc));
266 pc->data = (void *)lsc;
267
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 }
278